Detecting entanglement in quantum many-body systems via permutation moments

Multipartite entanglement plays an essential role in both quantum information science and many-body physics. Due to the exponentially large dimension and complex geometric structure of the state space, the detection of entanglement in many-body systems is extremely challenging in reality. Conventional means, like entanglement witness and entropy criterion, either highly depend on the prior knowledge of the studied systems or the detection capability is relatively weak. In this work, we propose a framework for designing multipartite entanglement criteria based on permutation moments, which have an effective implementation with either the generalized control-SWAP quantum circuits or the random unitary techniques. As an example, in the bipartite scenario, we develop an entanglement criterion that can detect bound entanglement and show strong detection capability in the multi-qubit Ising model with a long-range $XY$ Hamiltonian. In the multipartite case, the permutation-moment-based criteria can detect entangled states that are not detectable by any criteria extended from the bipartite case. Our framework also shows potential in entanglement quantification and entanglement structure detection.

Zhenhuan Liu

Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, China

Yifan Tang
Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, China
Department of Mathematics and Computer Science, Freie Universität Berlin, 14195 Berlin, Germany
Department of Physics, Freie Universität Berlin, 14195 Berlin, Germany

Hao Dai
Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, China

Pengyu Liu
Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, China

Shu Chen
Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, China

Xiongfeng Ma
Center for Quantum Information, Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing 100084, China

The past decades have witnessed great progress in understanding quantum entanglement [Horodecki et al., 2009]. To date, entanglement acts not only as the cornerstone of quantum information science, but also as a new perspective in many other fields, like quantum thermodynamics [Ueda, 2020, Kaufman et al., 2016], condensed matter physics [Laflorencie, 2016], and quantum gravity [Nishioka, 2018]. Especially in many-body physics, the dynamical behavior, scaling property, and spectral form of entanglement are key indicators to characterize different phases of the system [Abanin et al., 2019, Amico et al., 2008].

As a resource that cannot be produced by local operations and classical communication (LOCC), entangled $k$-partite states are those that cannot be written in a separable form, $\rho=\sum_ip_i\rho_{1}^i\otimes \cdots \otimes\rho_{k}^i$, where $p_i\ge 0$ satisfying the normalization condition $\sum_ip_i=1$ and $\rho_{r}^i$ is the density matrix of the $r$-th subsystem. As the dimension of a quantum system grows exponentially with the number of qubits, the geometric structure of state space becomes highly complicated, making entanglement detection a resource-consuming task. In fact, determining whether a state is entangled or not is generally a NP-hard problem [Gurvits, 2003]. For pure states or states with enough prior knowledge, entanglement can be effectively detected by purity measurements [Kaufman et al., 2016, Brydges et al., 2019], variational algorithms [Kokail et al., 2021], or entanglement witness [Gühne and Tóth, 2009, Zhang et al., 2021]. While in the noisy intermediate-scale quantum era [Preskill, 2018], the processed states are usually disturbed by unpredictable noise, rendering the detection capability of conventional means ineffective. Consequently, it is important to find implementable and efficient methods to detect multipartite entanglement with state-of-the-art devices.

For a generic mixed multipartite state without prior knowledge, there are two commonly-used techniques to detect entanglement, density matrix moments and index permutation. Moments of density matrix, $\tr(\rho^n)$, carry much information about the states and are relatively easy to measure [Ekert et al., 2002, van Enk and Beenakker, 2012, Elben et al., 2019].  Hence, they become practical tools in estimating properties of quantum systems [Smith et al., 2017], including quantum entanglement [Horodecki, 2003, Gühne and Tóth, 2009, Imai et al., 2021, Beckey et al., 2021, Ketterer et al., 2022]. However, most moment-based entanglement criteria are specially designed for states with few parties or low dimensions.
A general moment-based entanglement-detection framework for multipartite systems is still missing.

Based on the rearrangement of density matrix elements, the index permutation criterion [Horodecki et al., 2006] can be applied in systems with an arbitrary number of parties and dimensions and has many generalizations [Zhang et al., 2008]. In general, a $k$-partite quantum state can be represented using a matrix with $2k$ indices,

\begin{equation}
\begin{split}
\rho=\sum_{s_1,\cdots,s_{2k}}\rho_{s_1s_2,\cdots , s_{2k-1}s_{2k}}\ketbra{s_1\cdots s_{2k-1}}{s_{2}\cdots s_{2k}},
\end{split}
\end{equation}

where $s_1,s_3,\dots,s_{2k-1}$ represent the row indices, and $s_2,s_4,\dots,s_{2k}$ represent the column ones. The two indices, $s_{2r-1}$ and $s_{2r}$, denote for the $r$-th subsystem. By changing the order of these $2k$ indices, one gets a new matrix, $\mathcal{R}_{\pi}(\rho)$, with

\begin{equation}\label{eq:defRpi}
\begin{split}
\left[\mathcal{R}_{\pi}(\rho)\right]_{s_1s_2,\cdots, s_{2k-1}s_{2k}} = \rho_{s_{\pi(1)}s_{\pi(2)},\cdots , s_{\pi(2k-1)}s_{\pi(2k)}},
\end{split}
\end{equation}

where $\pi$ is an element of $2k$-th permutation group $\mathcal{S}_{2k}$. For simplicity, hereafter, we use $\mathcal{R}_\pi$ to denote $\mathcal{R}_\pi(\rho)$. Using the property of index permutation, one could prove that [Horodecki et al., 2006]

\begin{equation}\label{eq:permutationcriterion}
\begin{split}
\norm{\mathcal{R}_\pi}=\tr\left(\sqrt{\mathcal{R}_\pi\mathcal{R}_\pi^\dagger}\right)=\sum_i\lambda_i \le 1,
\end{split}
\end{equation}

for all $k$-partite separable states, where $\{\lambda_i\}$ are the singular values of $\mathcal{R}_\pi$. The violation of this inequality indicates entanglement.  In the bipartite scenario,  when setting $\pi$ to be $(1,2)$ and $(2,3)$, where $(\cdot,\cdot)$ denotes exchanging two indices, one gets the widely-used positive partial transposition (PPT) criterion [Peres, 1996] and the computable cross norm (CCNR) criterion [Chen and Wu, 2003], respectively. However, because index permutation is an unphysical operation, and the permutation criteria are based on singular value decomposition, a highly nonlinear operation, measurement of $\norm{\mathcal{R}_\pi}$ usually requires full state tomography, which is extremely resource-consuming [Haah et al., 2017].

To harness the power of permutation criteria in multipartite entanglement detection, we borrow the idea from moment criteria. Although it is generally hard to measure $\norm{\mathcal{R}_\pi}$ directly, one can alternatively estimate the higher-order moments, $M_{2n}^\pi=\tr[(\mathcal{R}_\pi\mathcal{R}_\pi^\dagger)^n]=\sum_i\lambda_i^{2n}$, which are much easier to access. These permutation moments can help to lower bound $\norm{\mathcal{R}_\pi}=\sum_i\lambda_i$ and infer whether the state is multipartite entangled or not. A similar idea has also been used in the estimation of quantum negativity [Gray et al., 2018, Yu et al., 2021] and entropy [Smith et al., 2017]. By changing the index permutation operation $\mathcal{R}_\pi(\cdot)$ and measuring different orders of moments, we generate a series of implementable multipartite entanglement criteria, which we call moment-based permutation criterion. The entanglement detection flowchart is shown in Fig. <ref>.

Flowchart of entanglement detection. To detect multipartite entanglement of $\rho$, one needs to first choose an index permutation operation $\mathcal{R}_\pi(\cdot)$ and set $\norm{\mathcal{R}_\pi(\rho)}=\sum_i\lambda_i$ as the entanglement indicator. Then, one measures the permutation moments $\{M_{2n}^\pi=\sum_i\lambda_i^{2n}\}_n$ and use these moments to lower bound $\norm{\mathcal{R}_\pi(\rho)}$. If the lower bound is larger than the entanglement threshold set for $\norm{\mathcal{R}_\pi(\rho)}$, the multipartite entanglement is successfully detected. Otherwise, one can measure higher-order moments or pick another index permutation and repeat the procedure.

Moment-based permutation criteria. —  For multipartite quantum state $\rho$, each party has two indices, one for row and one for column. Given $\mathcal{R}_\pi$, all the parties can be divided into four types, based on the position transition of their two indices, as shown in the first line of Table <ref>.

We find two properties of index permutation. First, if a party is T1-type or T2-type for $\mathcal{R}_\pi$, then it will keep the type for $\mathcal{R}_\pi^\dagger$; while if it is R1-type or R2-type for $\mathcal{R}_\pi$, then it becomes R2-type or R1-type for $\mathcal{R}_\pi^\dagger$, respectively. Second, the indices contraction in $M_{2n}^\pi=\tr[(\mathcal{R}_\pi\mathcal{R}_\pi^\dagger)^n]$ only acts on the indices from the same party of the $2n$ copies of $\rho$. Hence, if we list these $2n$ copies of states in order, the indices from an R-type party will contract with one of its two neighboring states and the indices from a T-type party will contract with both of its neighboring states. Note that $\mathcal{R}_\pi$ might not be a square matrix, so its odd moments are generally inaccessible.

Type     T1     T2     R1     R2

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho) at (0,0) $\rho$;
[Lline] (rho.west) –++ (-.2,0);
[Rline] (rho.east) –++ (.2,0);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho) at (0,0) $\rho$;
[Lline] (rho.west) –++ (-.15,0) to [bend left=80] ($(rho.east)+(.225,0)$) –++ (.15,0);
[Rline] (rho.east) –++ (.15,0) to [bend left=80] ($(rho.west)+(-.225,0)$) –++ (-.15,0);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho) at (0,0) $\rho$;
[Lline] (rho.west) –++ (-.3,0);
[Rline] (rho.east) –++ (.15,0) to [bend left=50] ($(rho.west)+(-.15,-0.15)$) –++ (-.15,0);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho) at (0,0) $\rho$;
[Lline] (rho.west) –++ (-.15,0) to [bend left=50] ($(rho.east)+(.15,0.15)$) –++ (.15,0);;
[Rline] (rho.east) –++ (.3,0);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.175) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0) node (L1) ;
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0) node (L2) ;
[Rline] (rho2.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.175) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0) node (L1) ;
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0) node (L2) ;
[Rline] (rho2.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.175) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0) node (L1) ;
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0) node (L2) ;
[Rline] (rho2.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.175) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0) node (L1) ;
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0) node (L2) ;
[Rline] (rho2.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.525) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0);
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0);
[Rline] (rho2.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho3) at ($(rho2)+(0,-.35)$) $\rho$;
[Lline] (rho3.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0) node (L1) ;
[Rline] (rho3.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho4) at ($(rho3)+(0,-.35)$) $\rho$;
[Lline] (rho4.west) –++ (-.1,0) –++ (-.2,1.05) –++ (-0.1,0) node (L2) ;
[Rline] (rho4.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.525) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-1.05) –++ (-0.1,0) node (L1) ;
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0) node (L2) ;
[Rline] (rho2.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho3) at ($(rho2)+(0,-.35)$) $\rho$;
[Lline] (rho3.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0);
[Rline] (rho3.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho4) at ($(rho3)+(0,-.35)$) $\rho$;
[Lline] (rho4.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0);
[Rline] (rho4.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.525) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-1.05) –++ (-0.1,0) node (L1) ;
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0);
[Rline] (rho2.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho3) at ($(rho2)+(0,-.35)$) $\rho$;
[Lline] (rho3.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0);
[Rline] (rho3.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho4) at ($(rho3)+(0,-.35)$) $\rho$;
[Lline] (rho4.west) –++ (-.1,0) –++ (-.2,1.05) –++ (-0.1,0) node (L2) ;
[Rline] (rho4.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

Lline/.style=-,blue,very thick,rounded corners,
Rline/.style=-,red,very thick,rounded corners,
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho1) at (0,.525) $\rho$;
[Lline] (rho1.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0);
[Rline] (rho1.east) –++ (.1,0) node (R1) ;
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho2) at ($(rho1)+(0,-.35)$) $\rho$;
[Lline] (rho2.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0) node (L1) ;
[Rline] (rho2.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho3) at ($(rho2)+(0,-.35)$) $\rho$;
[Lline] (rho3.west) –++ (-.1,0) –++ (-.2,-.35) –++ (-0.1,0) node (L2) ;
[Rline] (rho3.east) –++ (.1,0);
[draw=black,rectangle,rounded corners=0.1cm,minimum size=.42cm,,text=black,very thick] (rho4) at ($(rho3)+(0,-.35)$) $\rho$;
[Lline] (rho4.west) –++ (-.1,0) –++ (-.2,.35) –++ (-0.1,0);
[Rline] (rho4.east) –++ (.1,0) node (R2) ;
[Bline] (L1.south) – (L2.north);
[Bline] (R1.north) – (R2.south);

We use the tensor network to illustrate the four kinds of parties and their second and fourth moments. The boxes represent the subsystems of a generic $k$-partite state, and the two legs represent the row and column indices. The grey dashed lines represent periodic boundary condition. From the second line, one can find that the operators to estimate the second moments for all these four kinds of parties are the SWAP operators, which are represented by changing the order of two legs. Hence, $\tr(\mathcal{R}_\pi\mathcal{R}_\pi^\dagger)=\tr(\rho^2)$ for all $\pi\in\mathcal{S}_{2k}$. From the third line, one can find that the operators to estimate $\tr[(\mathcal{R}_\pi\mathcal{R}_\pi^\dagger)^2]$ for R-type parties are still SWAP operators, while the operators for T-type parties are cyclic permutation operators, which are represented by changing the order of the four legs cyclically.

Based on these results, we can prove that $M_{2n}^\pi$ can be measured directly by a joint observable.

Given a $k$-partite state $\rho$ and the index permutation operation $\mathcal{R}_\pi$, the $2n$-th moment of $\mathcal{R}_{\pi}$, $M_{2n}^\pi:=\tr[(\mathcal{R}_{\pi}\mathcal{R}_{\pi}^\dagger)^n]$, can be estimated by observable measurement on 2n copies of $\rho$,

\begin{equation}\label{eq:main res}
\begin{split}
M_{2n}^\pi=\tr\left(O^\pi_{2n}\rho^{\otimes 2n}\right)=\frac{1}{2}\tr\left[\left(\bigotimes_{i=1}^k U^\pi_{i}+h.c.\right) \rho^{\otimes 2n}\right].
\end{split}
\end{equation}

For T1-type parties $U_{i}^\pi=\overrightarrow{\Pi}_i$ and for T2-type parties $U_{i}^\pi=\overleftarrow{\Pi}_i$.
Here $\overrightarrow{\Pi}$ and $\overleftarrow{\Pi}$ are the cyclic permutation operators in different directions, satisfying $\overrightarrow{\Pi}\ket{s_1,\cdots,s_{2n}}=\ket{s_{2n},s_1,\cdots,s_{2n-1}}$ and $\overleftarrow{\Pi}\ket{s_1,\cdots,s_{2n}}=\ket{s_2,\cdots,s_{2n},s_1}$. For R1-type parties $
U_{i}^\pi=\mbb{S}^{(2n,1)}_i\otimes \mbb{S}^{(2,3)}_i\otimes \cdots \otimes \mbb{S}^{(2n-2,2n-1)}_i$ and for R2-type parties $U_{i}^\pi=\mbb{S}^{(1,2)}_i\otimes \mbb{S}^{(3,4)}_i\otimes \cdots \otimes \mbb{S}^{(2n-1,2n)}_i$, where $\mbb{S}^{(u,v)}$ is the SWAP operator acting on the $u$-th and $v$-th copies.

We leave the proofs of theorems in Appendix <ref>. The special cases of partial transposed and realigned moments for a two-qubit system have been discussed in Ref. [Carteret, 2005, Cai and Song, 2008].

Borrowing the ideas from [Ekert et al., 2002, Knill and Laflamme, 1998, Cai and Song, 2008], by introducing an ancilla qubit, we can design a quantum circuit to measure $M_{2n}^\pi$ based on the control-unitary operations, see Fig. <ref>. As the SWAP operators are the generators of the permutation group, all the control-unitary operators in this circuit can be decomposed into a polynomial number of the 3-qubit control-SWAP operators.

Quantum circuit for measuring $M_{2n}^\pi$. Inputs of this algorithm are $2n$ identical copies of $\rho$ and an ancilla qubit. The quantum gates in this circuit include Hadamard gates, labeled with $H$, and the control-$U_i^\pi$ gate. As $U_i^\pi$ can be decomposed to a series of SWAP operators, denoted by the connected $X$-form signs, this circuit can be realized by 3-qubit quantum gates. The measurement of the ancilla qubit is in computational basis.

However, the simultaneous preparation of $2n$ identical copies of $\rho$ is greatly challenging for state-of-the-art quantum devices. Fortunately, the recently developed techniques, shadow estimation [Aaronson, 2019, Huang et al., 2020] and randomized measurements [van Enk and Beenakker, 2012, Brydges et al., 2019, Elben et al., 2019], provide means to measure these moments by local (single-qubit) or global (multi-qubit) single-copy operations. Practically speaking, global operations are still challenging, and local protocols are the ones commonly used in real experiments. Shadow estimation has a wide range of applications while inefficient in general. Randomized measurements has lower sample complexities, while it can only measure some specialized physical quantities. In Appendix <ref>, we propose the measurement of the permutation moments using these two protocols and a hybrid one and analyze their sample complexities.

In the bipartite scenario, $M_n^{(1,2)}=\tr[(\mathcal{R}_{(1,2)})^n]=\tr[(\rho_{AB}^{T_A})^n]$ and $M_{2n}^{(2,3)}$ are key quantities that help to construct the weak-form PPT criteria [Elben et al., 2020, Yu et al., 2021, Neven et al., 2021] and the criteria proposed later in Eq. (<ref>) and Eq. (<ref>), respectively. The local randomized measurements scheme is not applicable for the measurement of $M_n^{(1,2)}$ [Zhou et al., 2020]. The existing single-copy local protocol for measuring $M_n^{(1,2)}$ requires the shadow scheme. While, according to Theorem <ref>, the observables for measuring $M_{2n}^{(2,3)}$ have a simple form. Thus, $M_{2n}^{(2,3)}$ can be measured through the local randomized measurements protocol and have a much lower sample complexity. We list the sample complexities of measuring $M_3^{(1,2)}$ and $M_{4}^{(2,3)}$ in Table <ref>.

Global Protocol     Local Protocol

$M_3^{(1,2)} $     $O(D^{\frac{2}{3}})$[Zhou et al., 2020]     $O(D^2)$[Elben et al., 2020]

$M_4^{(2,3)}$     $O(D^{\frac{1}{2}})$     $O(D^{1.187})$

This table shows the best-known sample complexities of measuring $M_3^{(1,2)}$ and  the complexities of protocols developed in Appendix <ref> to measure $M_4^{(2,3)}$. $D$ is the dimension of the underlying Hilbert space. This table shows a nearly quadratic improvement in the local case.

To find the lower bound of $\norm{\mathcal{R}_\pi}=\sum_i\lambda_i$ using these moments, one needs to solve an optimization problem. The original optimization problem is extremely hard to solve because we have an exponentially large number of $\lambda_i$.  Adopting the Lagrange multiplier method, we can simplify the optimization to a problem of solving a set of polynomial equations.

The minimum value of $\norm{\mathcal{R}_\pi}$ given $M_2^\pi$, ..., $M_{2n}^\pi$ is reached when there are at most $n$ non-zero $\lambda_i$s. Thus, denote the solution of this problem to be $E_{2n}^\pi(\rho)$, it is equal to the solution of the following optimization problem,

\begin{equation}
\begin{split}
\min_{q_1,\cdots,q_n\in\mathbb{N}}& E_{2n}^\pi(\rho)= q_1\lambda_1+q_2\lambda_2+\cdots+q_n\lambda_n\\
\mathrm{s.t.}& \ \sum_{i=1}^nq_i\lambda_i^2=M_2^\pi,\cdots, \ \sum_{i=1}^nq_i\lambda_i^{2n}=M_{2n}^\pi\\
& \ q_1+q_2+\cdots+q_n\le L,
\end{split}
\end{equation}

where $L$ is the number of the singular values of $\mathcal{R}_\pi$ and $q_i$ is the degeneracy of singular value $\lambda_i$.

As a special case, when we only know the value of $M_2^\pi$ and $M_4^\pi$, the minimum of $\sum_i\lambda_i$ has an analytical form [Berry and Sanders, 2003]

\begin{equation}\label{eq:E_4 definition}
\begin{split}
\end{split}
\end{equation}

where $q=\lfloor \frac{(M_2^\pi)^2}{M_4^\pi}\rfloor$ and $U=\sqrt{q(q+1)M_4^\pi-q (M_2^\pi)^2}$.

Now, we can formally represent the moment-based permutation criteria as

\begin{equation}\label{eq:criterionformulation}
\begin{split}
E_{2n}^\pi(\rho)\le 1 \ , \ \forall \pi\in\mathcal{S}_{2k} \ , \ n\in\mathbb{N}
\end{split}
\end{equation}

for all separable $k$-partite state $\rho$. In fact, $E_{2n}^\pi(\cdot)$ may not necessarily be the function of $\rho$. Adopting the bipartite entanglement criterion introduced in Ref. [Zhang et al., 2008], we get

\begin{equation}\label{eq:enhanced CCNR}
\begin{split}
E_{2n}^{(2,3)}(\rho_{AB}-\rho_A\otimes\rho_B)\le \sqrt{(1-\tr\rho_A^2)(1-\tr\rho_B^2)}
\end{split}
\end{equation}

for separable $\rho_{AB}$.

Bipartite entanglement detection. — Compared with existing entanglement detection schemes based on partial transposed moments [Elben et al., 2020, Yu et al., 2021, Neven et al., 2021], this framework is not only a direct generalization to multipartite entanglement, but also enhances the detection capability in the bipartite scenario.  With the second and fourth moments only, Eq. (<ref>) can detect $3\times 3$-dimensional bound entanglement constructed using the unextendible product basis proposed in Ref. [Bennett et al., 1999]. We leave the detailed discussion in Appendix <ref>.

The criterion of Eq. (<ref>) also performs well in practical physical systems. We study the local bipartite entanglement dynamics in a quantum system evolved under a long-range $XY$ Hamiltonian. Specifically, we choose a 10-qubit open boundary Ising model with the Hamiltonian of the form

\begin{equation}\label{eq:XY Hamiltonian}
\begin{split}
H_{XY}=\sum_{i<j}J_{ij}(\hat{\sigma}_i^+\hat{\sigma}_j^-+\hat{\sigma}_i^-\hat{\sigma}_j^+)+ B_z\sum_i\hat{\sigma}_i^z,
\end{split}
\end{equation}

where $\hat{\sigma}_i^z$, $\hat{\sigma}_i^+$, and $\hat{\sigma}_i^-$ are the spin-$\frac{1}{2}$ Pauli-$Z$, raising, and lowering operator acting on the $i$-th qubit; $J_{ij}=\frac{J_0}{|i-j|^\alpha}$ is the interaction strength following the power-law decay with $J_0$ and $\alpha$ set to be $420s^{-1}$ and $1.24$, respectively [Brydges et al., 2019]; $B_z$ stands for transverse field and is set to be $400s^{-1}$. This Hamiltonian has been realized in real physical systems [Jurcevic et al., 2014, Brydges et al., 2019] and is often served as the benchmark of detection capabilities of entanglement criteria [Neven et al., 2021, Elben et al., 2020].

Local entanglement decay in thermal system (Color Online). The entanglement dynamics of the local systems $A$ and $B$, which are marked by squares and initialized to be $\ket{\psi(t=0)}_{AB}=\frac{1}{\sqrt{2}}(\ket{0}^{\otimes N_{AB}}+\ket{1}^{\otimes N_{AB}})$. Qubits without squares are initialized to be the tensor product of $\ketbra{0}$, act as part $C$. The entanglement of $AB$ is detected when the value is above zero for each criterion. The grey areas represent the time periods in which the entanglement can only be detected by the $E_4^*$ criterion. }
\label{fig:Quench Dynamics}
\end{figure}

The 10-qubit chain is divided into three parts, $A$, $B$ and $C$, where $A$ and $B$ constitute the local system we study, initialized to be $1/√(2)(|0⟩^⊗N_AB+|1⟩^⊗N_AB)$. $C$ acts as the bath, which is initialized to be the tensor product of $|0⟩$. We compare four implementable nonlinear criteria in investigating the entanglement dynamics of systems composed of $A$ and $B$. The first two criteria are \comments{moment-based permutation criteria based on $M_4^{(2,3)}$, }  Eq.~\eqref{eq:criterionformulation} and Eq.~\eqref{eq:enhanced CCNR}, when setting $π=(2,3)$ and $n=2$, labeled by $E_4^(2,3)$ and $E_4^*$, respectively. Others are the entropy criterion based on the comparison of the purities [\citenamefont{Horodecki et~al.}, 2009], labeled by $P_2$; and the weak-form PPT criterion based on $M_3^(1,2)=[(ρ_AB^T_A)^3]$ [\citenamefont{Yu et~al.}, 2021, \citenamefont{Neven et~al.}, 2021], labeled by $P_3$. We define four quantities to represent these criteria which satisfy that $E(ρ)>0$ iff the entanglement is detected by the corresponding criterion.

The numerical simulation results [\citenamefont{Johansson et~al.}, 2012] are shown in Fig.~\ref{fig:Quench Dynamics}. One could find that \comments{the decay of the entanglement initially localized in $A$ and $B$ is observed by all four criteria. Besides, }the moment-based permutation criteria, especially $E_4^*$, have an obvious advantage since it detects entanglement while all others fail in various time periods and different choices of $A$ and $B$.

In addition to the strong detection capability, the key quantities in this framework, $E_2n^π(ρ)$, have clear mathematical meaning as they give the lower bounds of the permutation norms. The permutation norms, including entanglement negativity, can be treated as entanglement measures. We thus conjecture that these quantities can also be used as entanglement measures. In Appendix \ref{sec:physics}, we support this conjecture by showing that $E_4^(2,3)(ρ)$ can witness the entanglement scaling transition in a quantum dynamical phase transition [\citenamefont{Abanin et~al.}, 2019, \citenamefont{Serbyn et~al.}, 2013, \citenamefont{Smith et~al.}, 2016, \citenamefont{Wu and Das~Sarma}, 2016] and the entanglement rainbow structure for the eigenstates of a thermal Hamiltonian [\citenamefont{Ueda}, 2020, \citenamefont{Kim et~al.}, 2014, \citenamefont{Cotler et~al.}, 2021].

\emph{Multipartite entanglement detection.}---Another advantage of our framework lies in multipartite entanglement detection. There exist multipartite entangled states that are separable in any bipartition and thus cannot be detected by any criteria extended from the bipartite case [\citenamefont{Jungnitsch et~al.}, 2011], including the PPT and CCNR criteria.
Theorem~\ref{thm:LagrangeMultipler} provides us new means to design practical entanglement criteria for these states. We depict the sets of detectable multipartite entangled states of different criteria in Fig~\ref{fig:entlevel}.

An important example is also based on the unextendible product basis [\citenamefont{Bennett et~al.}, 1999]. Consider a three-qubit system and define four product pure states
\begin{equation}
\{\ket{\psi}_i\}_{i=1}^4=\{\ket{0,1,+},\ket{1,+,0},\ket{+,0,1},\ket{-,-,-}\},
\end{equation}
where $|±⟩=(|0⟩±|1⟩)/√(2)$. It has been proved that the state
\begin{equation}
\rho = \frac{1}{4}\left(\mathbb{I}_8-\sum_{i=1}^4\ketbra{\psi_i}{\psi_i}\right)
\end{equation}
is separable in any bipartition and thus its detection needs a new kind of moment-based permutation criterion other than PPT and CCNR. We find that when setting $π=1,2,3,4,5,61,3,2,4,5,6$, realigning the first two parties and keeping the third party unchanged, the entanglement of this state can be detected using $E_8^π(ρ)$, which only requires four orders of moments. We leave some details of calculating $E_8^π(ρ)$ in Appendix \ref{subsec:multientdetection}.

\begin{figure}
\centering
\includegraphics[width=0.45\textwidth]{EntanglementLevel.pdf}
\caption{Illustration of different sets of detectable states (Color Online). Bipartition PPT: PPT in any bipartition; Bipartition Separable: separable in any bipartition; Fully Separable: $\sum_ip_i\rho_{1}^i\otimes \cdots \otimes\rho_{k}^i$; Index Permutation: states that cannot be detected by an index permutation criterion other than the bipartite partial transposition; Moment Permutation: states that cannot be detected using finite numbers of permutation moments; Example State: a state that is separable in any bipartition while can be detected by a moment-based permutation criterion. When the number of moments increases, the dashed red circle will approach the solid red circle.
\label{fig:entlevel}
\end{figure}

For multipartite quantum systems, entanglement can have a rather complex entanglement structure [\citenamefont{Lu et~al.}, 2018]. At the same time, the tools for detecting entanglement structure are quite restrictive [\citenamefont{Ren et~al.}, 2021]. In Appendix \ref{sec:structure}, we show that our framework can also be generalized to detect the multipartite entanglement structure.

\emph{Outlook.}---The techniques we developed in this work, including the moment measurements and bounding the lower-order moment using the higher-order moments, have many potential applications, like the positive map entanglement detection [\citenamefont{Gühne and Tóth}, 2009] and the trace distance estimation. Furthermore, it is also interesting to investigate how to generalize the framework to entanglement detection in continuous variable systems [\citenamefont{Zhang et~al.}, 2013].

We thank Junjie Chen, Zhaohui Wei, and Xiaodong Yu for valuable discussions. This work was supported by the National Natural Science Foundation of China Grants No. 12174216 and No. 11875173 and the National Key Research and Development Program of China Grants No. 2019QY0702 and No. 2017YFA0303903.

\bibliographystyle{apsrev}
\begin{thebibliography}{57}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
\def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
\def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
\def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
\def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

[\citenamefont{Horodecki et~al.}, 2009]
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Horodecki}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Horodecki}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Horodecki}},
\bibnamefont{and}
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Horodecki}},
\bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{81}},
\bibinfo{pages}{865} (\bibinfo{year}{2009}),
\urlprefix\url{https://link.aps.org/doi/10.1103/RevModPhys.81.865}.

[\citenamefont{Ueda}, 2020]
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Ueda}},
\bibinfo{journal}{Nature Reviews Physics} \textbf{\bibinfo{volume}{2}},
\bibinfo{pages}{669} (\bibinfo{year}{2020}),
\urlprefix\url{https://www.nature.com/articles/s42254-020-0237-x.pdf}.

[\citenamefont{Kaufman et~al.}, 2016]
\bibinfo{author}{\bibfnamefont{A.~M.} \bibnamefont{Kaufman}},
\bibinfo{author}{\bibfnamefont{M.~E.} \bibnamefont{Tai}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Lukin}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Rispoli}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Schittko}},
\bibinfo{author}{\bibfnamefont{P.~M.} \bibnamefont{Preiss}},
\bibnamefont{and} \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Greiner}},
\bibinfo{journal}{Science} \textbf{\bibinfo{volume}{353}},
\bibinfo{pages}{794} (\bibinfo{year}{2016}),
\urlprefix\url{https://www.science.org/doi/abs/10.1126/science.aaf6725}.

[\citenamefont{Laflorencie}, 2016]
\bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Laflorencie}},
\bibinfo{journal}{Physics Reports} \textbf{\bibinfo{volume}{646}},
\bibinfo{pages}{1} (\bibinfo{year}{2016}), ISSN \bibinfo{issn}{0370-1573},
\bibinfo{note}{quantum entanglement in condensed matter systems},
\urlprefix\url{https://www.sciencedirect.com/science/article/pii/S0370157316301582}.

[\citenamefont{Nishioka}, 2018]
\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Nishioka}},
\bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{90}},
\bibinfo{pages}{035007} (\bibinfo{year}{2018}),
\urlprefix\url{https://link.aps.org/doi/10.1103/RevModPhys.90.035007}.

[\citenamefont{Abanin et~al.}, 2019]
\bibinfo{author}{\bibfnamefont{D.~A.} \bibnamefont{Abanin}},
\bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Altman}},
\bibinfo{author}{\bibfnamefont{I.}~\bibnamefont{Bloch}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Serbyn}},
\bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{91}},
\bibinfo{pages}{021001} (\bibinfo{year}{2019}),
\urlprefix\url{https://link.aps.org/doi/10.1103/RevModPhys.91.021001}.

[\citenamefont{Amico et~al.}, 2008]
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Amico}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Fazio}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Osterloh}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Vedral}},
\bibinfo{journal}{Rev. Mod. Phys.} \textbf{\bibinfo{volume}{80}},
\bibinfo{pages}{517} (\bibinfo{year}{2008}),
\urlprefix\url{https://link.aps.org/doi/10.1103/RevModPhys.80.517}.

[\citenamefont{Gurvits}, 2003]
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Gurvits}}, in
\emph{\bibinfo{booktitle}{Proceedings of the Thirty-Fifth Annual ACM
Symposium on Theory of Computing}} (\bibinfo{publisher}{Association for
Computing Machinery}, \bibinfo{address}{New York, NY, USA},
\bibinfo{year}{2003}), STOC '03, p. \bibinfo{pages}{10–19}, ISBN
\bibinfo{isbn}{1581136749},
\urlprefix\url{https://doi.org/10.1145/780542.780545}.

[\citenamefont{Brydges et~al.}, 2019]
\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Brydges}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Elben}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Jurcevic}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Vermersch}},
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Maier}},
\bibinfo{author}{\bibfnamefont{B.~P.} \bibnamefont{Lanyon}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zoller}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Blatt}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{C.~F.} \bibnamefont{Roos}},
\bibinfo{journal}{Science} \textbf{\bibinfo{volume}{364}},
\bibinfo{pages}{260} (\bibinfo{year}{2019}),
\urlprefix\url{https://www.science.org/doi/abs/10.1126/science.aau4963#pill-citations}.

[\citenamefont{Kokail et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Kokail}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Sundar}},
\bibinfo{author}{\bibfnamefont{T.~V.} \bibnamefont{Zache}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Elben}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Vermersch}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dalmonte}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{van Bijnen}},
\bibnamefont{and} \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zoller}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{127}},
\bibinfo{pages}{170501} (\bibinfo{year}{2021}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.127.170501}.

[\citenamefont{Gühne and Tóth}, 2009]
\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{Gühne}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Tóth}},
\bibinfo{journal}{Physics Reports} \textbf{\bibinfo{volume}{474}},
\bibinfo{pages}{1} (\bibinfo{year}{2009}), ISSN \bibinfo{issn}{0370-1573},
\urlprefix\url{https://www.sciencedirect.com/science/article/pii/S0370157309000623}.

[\citenamefont{Zhang et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Zhang}},
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Tang}},
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Zhou}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Ma}}, \bibinfo{journal}{Phys.
Rev. A} \textbf{\bibinfo{volume}{103}}, \bibinfo{pages}{052426}
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.103.052426}.

[\citenamefont{Preskill}, 2018]
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Preskill}},
\bibinfo{journal}{{Quantum}} \textbf{\bibinfo{volume}{2}},
\bibinfo{pages}{79} (\bibinfo{year}{2018}), ISSN \bibinfo{issn}{2521-327X},
\urlprefix\url{https://doi.org/10.22331/q-2018-08-06-79}.

[\citenamefont{Ekert et~al.}, 2002]
\bibinfo{author}{\bibfnamefont{A.~K.} \bibnamefont{Ekert}},
\bibinfo{author}{\bibfnamefont{C.~M.} \bibnamefont{Alves}},
\bibinfo{author}{\bibfnamefont{D.~K.~L.} \bibnamefont{Oi}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Horodecki}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Horodecki}},
\bibnamefont{and} \bibinfo{author}{\bibfnamefont{L.~C.} \bibnamefont{Kwek}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{88}},
\bibinfo{pages}{217901} (\bibinfo{year}{2002}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.88.217901}.

[\citenamefont{van Enk and Beenakker}, 2012]
\bibinfo{author}{\bibfnamefont{S.~J.} \bibnamefont{van Enk}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{C.~W.~J.} \bibnamefont{Beenakker}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{108}},
\bibinfo{pages}{110503} (\bibinfo{year}{2012}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.108.110503}.

[\citenamefont{Elben et~al.}, 2019]
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Elben}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Vermersch}},
\bibinfo{author}{\bibfnamefont{C.~F.} \bibnamefont{Roos}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zoller}},
\bibinfo{journal}{Phys. Rev. A} \textbf{\bibinfo{volume}{99}},
\bibinfo{pages}{052323} (\bibinfo{year}{2019}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.99.052323}.

[\citenamefont{Smith et~al.}, 2017]
\bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Smith}},
\bibinfo{author}{\bibfnamefont{J.~A.} \bibnamefont{Smolin}},
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Yuan}},
\bibinfo{author}{\bibfnamefont{Q.}~\bibnamefont{Zhao}},
\bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Girolami}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Ma}}, \bibinfo{journal}{arXiv
preprint arXiv:1707.09928}  (\bibinfo{year}{2017}),
\urlprefix\url{https://arxiv.org/pdf/1707.09928.pdf}.

[\citenamefont{Horodecki}, 2003]
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Horodecki}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{90}},
\bibinfo{pages}{167901} (\bibinfo{year}{2003}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.90.167901}.

[\citenamefont{Imai et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Imai}},
\bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Wyderka}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Ketterer}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{G\"uhne}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{126}},
\bibinfo{pages}{150501} (\bibinfo{year}{2021}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.126.150501}.

[\citenamefont{Beckey et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{J.~L.} \bibnamefont{Beckey}},
\bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Gigena}},
\bibinfo{author}{\bibfnamefont{P.~J.} \bibnamefont{Coles}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Cerezo}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{127}},
\bibinfo{pages}{140501} (\bibinfo{year}{2021}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.127.140501}.

[\citenamefont{Ketterer et~al.}, 2022]
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Ketterer}},
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Imai}},
\bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Wyderka}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{G\"uhne}},
\bibinfo{journal}{Phys. Rev. A} \textbf{\bibinfo{volume}{106}},
\bibinfo{pages}{L010402} (\bibinfo{year}{2022}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.106.L010402}.

[\citenamefont{Horodecki et~al.}, 2006]
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Horodecki}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Horodecki}},
\bibnamefont{and}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Horodecki}},
\bibinfo{journal}{Open Systems \& Information Dynamics}
\textbf{\bibinfo{volume}{13}}, \bibinfo{pages}{103} (\bibinfo{year}{2006}),
\urlprefix\url{https://link.springer.com/content/pdf/10.1007/s11080-006-7271-8.pdf}.

[\citenamefont{Zhang et~al.}, 2008]
\bibinfo{author}{\bibfnamefont{C.-J.} \bibnamefont{Zhang}},
\bibinfo{author}{\bibfnamefont{Y.-S.} \bibnamefont{Zhang}},
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Zhang}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{G.-C.} \bibnamefont{Guo}},
\bibinfo{journal}{Phys. Rev. A} \textbf{\bibinfo{volume}{77}},
\bibinfo{pages}{060301} (\bibinfo{year}{2008}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.77.060301}.

[\citenamefont{Peres}, 1996]
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Peres}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{77}},
\bibinfo{pages}{1413} (\bibinfo{year}{1996}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.77.1413}.

[\citenamefont{Chen and Wu}, 2003]
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Chen}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{L.-A.} \bibnamefont{Wu}},
\bibinfo{journal}{Quantum Info. Comput.} \textbf{\bibinfo{volume}{3}},
\bibinfo{pages}{193–202} (\bibinfo{year}{2003}), ISSN
\bibinfo{issn}{1533-7146}.

[\citenamefont{Haah et~al.}, 2017]
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Haah}},
\bibinfo{author}{\bibfnamefont{A.~W.} \bibnamefont{Harrow}},
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Ji}},
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Wu}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{N.}~\bibnamefont{Yu}}, \bibinfo{journal}{IEEE
Transactions on Information Theory} \textbf{\bibinfo{volume}{63}},
\bibinfo{pages}{5628} (\bibinfo{year}{2017}),
\urlprefix\url{https://ieeexplore.ieee.org/abstract/document/7956181}.

[\citenamefont{Gray et~al.}, 2018]
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Gray}},
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Banchi}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Bayat}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Bose}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{121}},
\bibinfo{pages}{150503} (\bibinfo{year}{2018}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.121.150503}.

[\citenamefont{Yu et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{X.-D.} \bibnamefont{Yu}},
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Imai}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{G\"uhne}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{127}},
\bibinfo{pages}{060504} (\bibinfo{year}{2021}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.127.060504}.

[\citenamefont{Carteret}, 2005]
\bibinfo{author}{\bibfnamefont{H.~A.} \bibnamefont{Carteret}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{94}},
\bibinfo{pages}{040502} (\bibinfo{year}{2005}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.94.040502}.

[\citenamefont{Cai and Song}, 2008]
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Cai}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Song}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{101}},
\bibinfo{pages}{190503} (\bibinfo{year}{2008}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.101.190503}.

[\citenamefont{Knill and Laflamme}, 1998]
\bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Knill}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Laflamme}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{81}},
\bibinfo{pages}{5672} (\bibinfo{year}{1998}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.81.5672}.

[\citenamefont{Aaronson}, 2019]
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Aaronson}},
\bibinfo{journal}{SIAM Journal on Computing} \textbf{\bibinfo{volume}{49}},
\bibinfo{pages}{STOC18} (\bibinfo{year}{2019}),
\urlprefix\url{https://epubs.siam.org/doi/pdf/10.1137/18M120275X}.

[\citenamefont{Huang et~al.}, 2020]
\bibinfo{author}{\bibfnamefont{H.-Y.} \bibnamefont{Huang}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Kueng}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Preskill}},
\bibinfo{journal}{Nature Physics} \textbf{\bibinfo{volume}{16}},
\bibinfo{pages}{1050} (\bibinfo{year}{2020}),
\urlprefix\url{https://doi.org/10.1038/s41567-020-0932-7}.

[\citenamefont{Elben et~al.}, 2020]
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Elben}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Kueng}},
\bibinfo{author}{\bibfnamefont{H.-Y.~R.} \bibnamefont{Huang}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{van Bijnen}},
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Kokail}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dalmonte}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Calabrese}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Kraus}},
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Preskill}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zoller}},
\bibnamefont{et~al.}, \bibinfo{journal}{Phys. Rev. Lett.}
\textbf{\bibinfo{volume}{125}}, \bibinfo{pages}{200501}
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.125.200501}.

[\citenamefont{Neven et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Neven}},
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Carrasco}},
\bibinfo{author}{\bibfnamefont{V.}~\bibnamefont{Vitale}},
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Kokail}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Elben}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Dalmonte}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Calabrese}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zoller}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Vermersch}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Kueng}},
\bibnamefont{et~al.}, \bibinfo{journal}{npj Quantum Information}
\textbf{\bibinfo{volume}{7}}, \bibinfo{pages}{152} (\bibinfo{year}{2021}),
ISSN \bibinfo{issn}{2056-6387},
\urlprefix\url{https://doi.org/10.1038/s41534-021-00487-y}.

[\citenamefont{Zhou et~al.}, 2020]
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Zhou}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zeng}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Liu}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{125}},
\bibinfo{pages}{200502} (\bibinfo{year}{2020}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.125.200502}.

[\citenamefont{Berry and Sanders}, 2003]
\bibinfo{author}{\bibfnamefont{D.~W.} \bibnamefont{Berry}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{B.~C.} \bibnamefont{Sanders}},
\bibinfo{journal}{Journal of Physics A: Mathematical and General}
\textbf{\bibinfo{volume}{36}}, \bibinfo{pages}{12255–12265}
(\bibinfo{year}{2003}), ISSN \bibinfo{issn}{1361-6447},
\urlprefix\url{http://dx.doi.org/10.1088/0305-4470/36/49/008}.

[\citenamefont{Bennett et~al.}, 1999]
\bibinfo{author}{\bibfnamefont{C.~H.} \bibnamefont{Bennett}},
\bibinfo{author}{\bibfnamefont{D.~P.} \bibnamefont{DiVincenzo}},
\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Mor}},
\bibinfo{author}{\bibfnamefont{P.~W.} \bibnamefont{Shor}},
\bibinfo{author}{\bibfnamefont{J.~A.} \bibnamefont{Smolin}},
\bibnamefont{and} \bibinfo{author}{\bibfnamefont{B.~M.}
\bibnamefont{Terhal}}, \bibinfo{journal}{Phys. Rev. Lett.}
\textbf{\bibinfo{volume}{82}}, \bibinfo{pages}{5385} (\bibinfo{year}{1999}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.82.5385}.

[\citenamefont{Jurcevic et~al.}, 2014]
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Jurcevic}},
\bibinfo{author}{\bibfnamefont{B.~P.} \bibnamefont{Lanyon}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Hauke}},
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Hempel}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zoller}},
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Blatt}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{C.~F.} \bibnamefont{Roos}},
\bibinfo{journal}{Nature} \textbf{\bibinfo{volume}{511}},
\bibinfo{pages}{202} (\bibinfo{year}{2014}),
\urlprefix\url{https://www.nature.com/articles/nature13461.pdf}.

[\citenamefont{Johansson et~al.}, 2012]
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Johansson}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Nation}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Nori}},
\bibinfo{journal}{Computer Physics Communications}
\textbf{\bibinfo{volume}{183}}, \bibinfo{pages}{1760} (\bibinfo{year}{2012}),
ISSN \bibinfo{issn}{0010-4655},
\urlprefix\url{https://www.sciencedirect.com/science/article/pii/S0010465512000835}.

[\citenamefont{Serbyn et~al.}, 2013]
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Serbyn}},
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Papi\ifmmode~\acute{c}\else
\'{c}\fi{}}}, \bibnamefont{and} \bibinfo{author}{\bibfnamefont{D.~A.}
\bibnamefont{Abanin}}, \bibinfo{journal}{Phys. Rev. Lett.}
\textbf{\bibinfo{volume}{111}}, \bibinfo{pages}{127201}
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.111.127201}.

[\citenamefont{Smith et~al.}, 2016]
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Smith}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Lee}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Richerme}},
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Neyenhuis}},
\bibinfo{author}{\bibfnamefont{P.~W.} \bibnamefont{Hess}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Hauke}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Heyl}},
\bibinfo{author}{\bibfnamefont{D.~A.} \bibnamefont{Huse}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Monroe}},
\bibinfo{journal}{Nature Physics} \textbf{\bibinfo{volume}{12}},
\bibinfo{pages}{907} (\bibinfo{year}{2016}),
\urlprefix\url{https://www.nature.com/articles/nphys3783.pdf}.

[\citenamefont{Wu and Das~Sarma}, 2016]
\bibinfo{author}{\bibfnamefont{Y.-L.} \bibnamefont{Wu}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Das~Sarma}},
\bibinfo{journal}{Phys. Rev. A} \textbf{\bibinfo{volume}{93}},
\bibinfo{pages}{022332} (\bibinfo{year}{2016}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.93.022332}.

[\citenamefont{Kim et~al.}, 2014]
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Kim}},
\bibinfo{author}{\bibfnamefont{T.~N.} \bibnamefont{Ikeda}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{D.~A.} \bibnamefont{Huse}},
\bibinfo{journal}{Phys. Rev. E} \textbf{\bibinfo{volume}{90}},
\bibinfo{pages}{052105} (\bibinfo{year}{2014}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevE.90.052105}.

[\citenamefont{Cotler et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{J.~S.} \bibnamefont{Cotler}},
\bibinfo{author}{\bibfnamefont{D.~K.} \bibnamefont{Mark}},
\bibinfo{author}{\bibfnamefont{H.-Y.} \bibnamefont{Huang}},
\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Hernandez}},
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Choi}},
\bibinfo{author}{\bibfnamefont{A.~L.} \bibnamefont{Shaw}},
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Endres}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Choi}},
\bibinfo{journal}{arXiv preprint arXiv:2103.03536}  (\bibinfo{year}{2021}),
\urlprefix\url{https://arxiv.org/pdf/2103.03536.pdf}.

[\citenamefont{Jungnitsch et~al.}, 2011]
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Jungnitsch}},
\bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Moroder}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{O.}~\bibnamefont{G\"uhne}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{106}},
\bibinfo{pages}{190502} (\bibinfo{year}{2011}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.106.190502}.

[\citenamefont{Lu et~al.}, 2018]
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Lu}},
\bibinfo{author}{\bibfnamefont{Q.}~\bibnamefont{Zhao}},
\bibinfo{author}{\bibfnamefont{Z.-D.} \bibnamefont{Li}},
\bibinfo{author}{\bibfnamefont{X.-F.} \bibnamefont{Yin}},
\bibinfo{author}{\bibfnamefont{X.}~\bibnamefont{Yuan}},
\bibinfo{author}{\bibfnamefont{J.-C.} \bibnamefont{Hung}},
\bibinfo{author}{\bibfnamefont{L.-K.} \bibnamefont{Chen}},
\bibinfo{author}{\bibfnamefont{L.}~\bibnamefont{Li}},
\bibinfo{author}{\bibfnamefont{N.-L.} \bibnamefont{Liu}},
\bibinfo{author}{\bibfnamefont{C.-Z.} \bibnamefont{Peng}},
\bibnamefont{et~al.}, \bibinfo{journal}{Phys. Rev. X}
\textbf{\bibinfo{volume}{8}}, \bibinfo{pages}{021072} (\bibinfo{year}{2018}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevX.8.021072}.

[\citenamefont{Ren et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Ren}},
\bibinfo{author}{\bibfnamefont{W.}~\bibnamefont{Li}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Smerzi}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Gessner}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{126}},
\bibinfo{pages}{080502} (\bibinfo{year}{2021}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.126.080502}.

[\citenamefont{Zhang et~al.}, 2013]
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Zhang}},
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Yu}},
\bibinfo{author}{\bibfnamefont{Q.}~\bibnamefont{Chen}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{C.~H.} \bibnamefont{Oh}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{111}},
\bibinfo{pages}{190501} (\bibinfo{year}{2013}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.111.190501}.

[\citenamefont{Wood et~al.}, 2011]
\bibinfo{author}{\bibfnamefont{C.~J.} \bibnamefont{Wood}},
\bibinfo{author}{\bibfnamefont{J.~D.} \bibnamefont{Biamonte}},
\bibnamefont{and} \bibinfo{author}{\bibfnamefont{D.~G.} \bibnamefont{Cory}},
\bibinfo{journal}{arXiv preprint arXiv:1111.6950}  (\bibinfo{year}{2011}),
\urlprefix\url{https://arxiv.org/pdf/1111.6950.pdf}.

[\citenamefont{Gu}, 2013]
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Gu}}, Ph.D. thesis
\urlprefix\url{https://qspace.library.queensu.ca/bitstream/handle/1974/8241/Gu_Yinzheng_201308_MSc.pdf?sequence=1}.

[\citenamefont{Roberts and Yoshida}, 2017]
\bibinfo{author}{\bibfnamefont{D.~A.} \bibnamefont{Roberts}} \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Yoshida}},
\bibinfo{journal}{Journal of High Energy Physics}
\textbf{\bibinfo{volume}{2017}}, \bibinfo{pages}{1} (\bibinfo{year}{2017}),
\urlprefix\url{https://link.springer.com/content/pdf/10.1007/JHEP04(2017)121.pdf}.

[\citenamefont{Zhu}, 2017]
\bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Zhu}}, \bibinfo{journal}{Phys.
Rev. A} \textbf{\bibinfo{volume}{96}}, \bibinfo{pages}{062336}
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.96.062336}.

[\citenamefont{Horn et~al.}, 1994]
\bibinfo{author}{\bibfnamefont{R.~A.} \bibnamefont{Horn}},
\bibinfo{author}{\bibfnamefont{R.~A.} \bibnamefont{Horn}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{C.~R.} \bibnamefont{Johnson}},
\emph{\bibinfo{title}{Topics in matrix analysis}}
(\bibinfo{publisher}{Cambridge university press}, \bibinfo{year}{1994}).

[\citenamefont{Li et~al.}, 2017]
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Li}},
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Wang}},
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Shen}},
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Chen}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{S.-M.} \bibnamefont{Fei}},
\bibinfo{journal}{Scientific Reports} \textbf{\bibinfo{volume}{7}},
\bibinfo{pages}{17274} (\bibinfo{year}{2017}), ISSN
\bibinfo{issn}{2045-2322},
\urlprefix\url{https://doi.org/10.1038/s41598-017-17585-7}.

[\citenamefont{Liu et~al.}, 2022]
\bibinfo{author}{\bibfnamefont{Z.}~\bibnamefont{Liu}},
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Zeng}},
\bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Zhou}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Gu}}, \bibinfo{journal}{Phys.
Rev. A} \textbf{\bibinfo{volume}{105}}, \bibinfo{pages}{022407}
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevA.105.022407}.

[\citenamefont{Rath et~al.}, 2021]
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Rath}},
\bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Branciard}},
\bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Minguzzi}}, \bibnamefont{and}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Vermersch}},
\bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{127}},
\bibinfo{pages}{260501} (\bibinfo{year}{2021}),
\urlprefix\url{https://link.aps.org/doi/10.1103/PhysRevLett.127.260501}.

\end{thebibliography}

\appendix

\onecolumngrid
\newpage

\setcounter{theorem}{0}
\setcounter{table}{0}

\section{Preliminaries}\label{sec:Preliminary}
\subsection{Tensor Network Basis}
As our work is mainly based on index permutation criterion, much tensor calculation is needed to derive the results. So here we will introduce a graphical method to conduct the tensor calculation, the tensor network [\citenamefont{Wood et~al.}, 2011], which also plays a crucial role in quantum simulation. In the following context, we will frequently use the technique introduced in this section to do the tensor calculation.

In tensor network representation, a matrix is represented as a box with open legs, which correspond to uncontracted indices. Those left-oriented legs stand for row indices and right-oriented ones stand for column indices. Vectors and scalars are represented as boxes with legs of same orientation and boxes with no leg, respectively. The number of legs depends on the number of parties we are interested in. A $k$-partite state is represented by a tensor with $k$ pairs of legs, with each pair of legs standing for row and column indices of each party. Connecting legs stands for index contraction, like matrix production and taking trace:
\begin{eqnarray}
\begin{tabular}{c}
\includegraphics[scale=0.15]{matrix_production.png}
\end{tabular}
, \ \ \
\tr(A)=
\begin{tabular}{c}
\includegraphics[scale=0.15]{traceA.png}
\end{tabular}.
\end{eqnarray}
For tensor production operation $A⊗B$, where no indices are contracted, $A$ and $B$ are just put together with no leg connection:
\begin{eqnarray}
A\otimes B=
\begin{tabular}{c}
\includegraphics[scale=0.15]{AtensorB.png}
\end{tabular}
\end{eqnarray}

Index permutation operations can be easily represented by changing the order of legs. Take partial transposition map and realignment map as examples,
\begin{eqnarray}
\rho_{AB}=
\begin{tabular}{c}
\includegraphics[scale=0.15]{rhoAB.png}
\end{tabular}
, \ \ \
\mathcal{R}_{(2,3)}(\rho_{AB})=
\begin{tabular}{c}
\includegraphics[scale=0.15]{R_rhoAB_.png}
\end{tabular}
, \ \ \
\rho_{AB}^{T_A}=\mathcal{R}_{(1,2)}(\rho_{AB})=
\begin{tabular}{c}
\includegraphics[scale=0.15]{rhoABTA.png}
\end{tabular}
\end{eqnarray}

In tensor network, permutation operators can be represented by changing the position of legs. SWAP operator is represented by exchanging two legs, cyclic operator is represented by sequentially moving each leg to the position of its neighboring leg. The operators we use to calculate the fourth moment of $ℛ_π$ can be graphically represented as:
\begin{eqnarray}
\mathbb{S}^{(1,2)}\otimes\mathbb{S}^{(3,4)}=
\begin{tabular}{c}
\includegraphics[scale=0.15]{SWAP1.png}
\end{tabular}
, \ \
\mathbb{S}^{(2,3)}\otimes\mathbb{S}^{(4,1)}=
\begin{tabular}{c}
\includegraphics[scale=0.15]{SWAP2.png}
\end{tabular}
, \ \
\overrightarrow{\Pi}=
\begin{tabular}{c}
\includegraphics[scale=0.15]{cyclic2.png}
\end{tabular}
, \ \
\overleftarrow{\Pi}=
\begin{tabular}{c}
\includegraphics[scale=0.15]{cyclic1.png}
\end{tabular}
\end{eqnarray}

\subsection{Random Unitary Basis}\label{subsec:random}
According to Schur-Weyl duality, random unitary is closely related to permutation operator, and is the basis of shadow estimation and randomized measurements, which will be discussed later. In unitary group, there exists a unique measure, which is called Haar measure, that satisfies
\begin{equation}
\begin{split}
\int_{\mathrm{Haar}}dU=1,\ \int_{\mathrm{Haar}}dUf(U)=\int_{\mathrm{Haar}}dUf(VU)=\int_{\mathrm{Haar}}dUf(UV)
\end{split}
\end{equation}
for arbitrary unitary $V$, where $U$ is integrated by Haar measure. Using Haar measure, one can define the $t$-fold twirling channel
\begin{equation}
\begin{split}
\Phi_t(O)=\int_{\mathrm{Haar}}dUU^{\otimes t}OU^{\dagger\otimes t},
\end{split}
\end{equation}
which equals to the linear combination of permutation operators [\citenamefont{Gu}, 2013, \citenamefont{Roberts and Yoshida}, 2017]
\begin{equation}
\begin{split}
\Phi_t(O)=\sum_{\pi,\sigma\in\mathcal{S}_t}C_{\pi,\sigma}\tr\left(\hat{W}_\pi O\right)\hat{W}_\sigma,
\end{split}
\end{equation}
where $𝒮_t$ is the permutation group of order $t$, $π$ and $σ$ are its elements, $C_π,σ$ is the element of Weigarten matrix, and $Ŵ_π$ is the permutation operator corresponding to $π$.

However, to construct a $t$-fold twirling channel, Haar measure may not be necessary, one can only average over a unitary ensemble $ℰ_t$ with a finite number of elements
\begin{equation}
\begin{split}
\Phi_t(O)=\frac{1}{|\mathcal{E}_t|}\sum_{U\in\mathcal{E}_t}U^{\otimes t}O U^{\dagger \otimes t},
\end{split}
\end{equation}
where $|ℰ_t|$ denotes the number of elements in $ℰ_t$. All the unitary ensembles that satisfy this equation are called unitary $t$-design. The commonly-used Clifford group has been proved to be a unitary 3-design for multi-qubit systems [\citenamefont{Zhu}, 2017].

\section{Proof of Theorems}\label{sec:proofs}

\subsection{Observable of Permutation Moments}\label{subsec:proofofobservable}

\begin{theorem}
Given a $k$-partite state $\rho$ and the index permutation operation $\mathcal{R}_\pi$, the $2n$-th moment of $\mathcal{R}_{\pi}$, $M_{2n}^\pi:=\tr[(\mathcal{R}_{\pi}\mathcal{R}_{\pi}^\dagger)^n]$, can be estimated by observable measurement on 2n copies of $\rho$,
\begin{equation}\label{eq:main res}
\begin{split}
M_{2n}^\pi=\tr\left(O^\pi_{2n}\rho^{\otimes 2n}\right)=\frac{1}{2}\tr\left[\left(\bigotimes_{i=1}^k U^\pi_{i}+h.c.\right) \rho^{\otimes 2n}\right].
\end{split}
\end{equation}
For T1-type parties $U_{i}^\pi=\overrightarrow{\Pi}_i$ and for T2-type parties $U_{i}^\pi=\overleftarrow{\Pi}_i$.
Here $\overrightarrow{\Pi}$ and $\overleftarrow{\Pi}$ are the cyclic permutation operators in different direction, satisfying $\overrightarrow{\Pi}\ket{s_1,\cdots,s_{2n}}=\ket{s_{2n},s_1,\cdots,s_{2n-1}}$ and $\overleftarrow{\Pi}\ket{s_1,\cdots,s_{2n}}=\ket{s_2,\cdots,s_{2n},s_1}$. For R1-type parties $
U_{i}^\pi=\mbb{S}^{(2n,1)}_i\otimes \mbb{S}^{(2,3)}_i\otimes \cdots \otimes \mbb{S}^{(2n-2,2n-1)}_i$ and for R2-type parties $U_{i}^\pi=\mbb{S}^{(1,2)}_i\otimes \mbb{S}^{(3,4)}_i\otimes \cdots \otimes \mbb{S}^{(2n-1,2n)}_i$, where $\mbb{S}^{(u,v)}$ is the SWAP operator acting on the $u$-th and $v$-th copies.
\end{theorem}

\begin{proof}
The key of this proof is to figure out the index permutation rules in $\tr[(\mathcal{R}_\pi\mathcal{R}_\pi^\dagger)^n]$. First recall the definition of index permutation operation:
\begin{equation}\label{eq:A1}
\begin{split}
\left[\mathcal{R}_\pi\right]_{s_1s_2,\cdots,s_{2k-1}s_{2k}}=\rho_{s_{\pi(1)}s_{\pi(2)},\cdots,s_{\pi(2k-1)}s_{\pi(2k)}}.
\end{split}
\end{equation}
According to the linearity of $\mathcal{R}_\pi(\cdot)$ and hermitian of $\rho$,
\begin{equation}
\begin{split}
\mathcal{R}_\pi^\dagger=[\mathcal{R}_\pi(\rho)^*]^T=\mathcal{R}_\pi(\rho^*)^T=\mathcal{R}_\pi(\rho^T)^T,
\end{split}
\end{equation}
where $T$ denotes the transposition operation. Hence, the element of $\mathcal{R}_\pi(\rho)^\dagger$ is
\begin{equation}\label{eq:A3}
\begin{split}
\left[\mathcal{R}_\pi^\dagger\right]_{s_1s_2,\cdots,s_{2k-1}s_{2k}}&=[\mathcal{R}_\pi(\rho^T)]_{s_2s_1,\cdots,s_{2k}s_{2k-1}}\\
\end{split}
\end{equation}
where $\pi'(\cdot)$ denotes the corresponding row number or column number of $\pi(\cdot)$, if $\pi(\cdot)$ is odd, then $\pi'(\cdot)=\pi(\cdot)+1$; otherwise, $\pi'(\cdot)=\pi(\cdot)-1$. Make a comparison of Eq.~\eqref{eq:A1} and Eq.~\eqref{eq:A3}, we find that the $r$-th row index of $\mathcal{R}_\pi$ and the $r$-th column index of $\mathcal{R}_\pi^\dagger$ are the row (column) and column (row) indices from the same subsystem of $\rho$. And this is same for the $r$-th column index of $\mathcal{R}_\pi$ and $r$-th row index of $\mathcal{R}_\pi^\dagger$. Using tensor network representation, this conclusion tells us that if transversely drawing the tensors of $\mathcal{R}_\pi$ and $\mathcal{R}_\pi^\dagger$, we will find they are mirror symmetric. We take $\mathcal{R}_{(2,3)}$, $\mathcal{R}_{(2,3)}$ and a generic permutation matrix $\mathcal{R}_\pi$ as examples:
\begin{eqnarray}
\mathcal{R}_{(2,3)}^\dagger=&
\left(
\begin{tabular}{c}
\includegraphics[scale=0.12]{R_rhoAB_.png}
\end{tabular}
\right)^\dagger=
\begin{tabular}{c}
\includegraphics[scale=0.12]{R_rhoAB_dagger.png}
\end{tabular}
, \ \ \
\mathcal{R}_{(1,2)}=
\left(
\begin{tabular}{c}
\includegraphics[scale=0.15]{rhoABTA.png}
\end{tabular}
\right)^\dagger=
\begin{tabular}{c}
\includegraphics[scale=0.15]{rhoABTAdagger.png}
\end{tabular}
\left(
\begin{tabular}{c}
\includegraphics[scale=0.12]{Rpirho.png}
\end{tabular}
\right)^\dagger=
\begin{tabular}{c}
\includegraphics[scale=0.12]{Rpirhodagger.png}
\end{tabular}.
\end{eqnarray}

Because of this mirror symmetry property, when multiplying $\mathcal{R}_\pi$ and $\mathcal{R}_\pi^\dagger$, leg connections only occur in the same subsystem. Besides, for those R1-type subsystems, their two legs connect in the multiplication $\mathcal{R}_\pi\mathcal{R}_\pi^\dagger$, and for those R2-type subsystems, their two legs connect in $\mathcal{R}_\pi^\dagger\mathcal{R}_\pi$. While for T-type subsystem, because its two legs are on both sides of $\mathcal{R}_\pi$ and $\mathcal{R}_\pi^\dagger$, so it will connect with both of its neighboring $\rho$. Hence, in $\tr[(\mathcal{R}_\pi\mathcal{R}_\pi^\dagger)^n]$, one would find that the R-type parties of two neighboring $\rho$ are either connected by legs or have no connection at all, depending on their exact types; all the T-type parties are connected by their legs, while the directions of connection are different for T1-type and T2-type parties. This connection rule is equivalent to measuring SWAP operators, $\mbb{S}$, on connected neighboring R-type subsystems and cyclic permutation operators $\overrightarrow{\Pi}$ or $\overleftarrow{\Pi}$ on $2n$ copies of T-type subsystems.

To summarize,
\begin{equation}\label{eq:M2nunitary}
\begin{split}
M_{2n}^\pi=\tr\left[\left(\bigotimes_{i=1}^kU_i^\pi\right)\rho^{\otimes 2n}\right]
\end{split}
\end{equation}
where $U^\pi_i=\mathbb{S}_i^{(2,3)}\otimes \cdots\otimes\mathbb{S}_i^{(2n,1)}$ for R1-type parties, $U^\pi_i=\mathbb{S}_i^{(1,2)}\otimes \cdots\otimes\mathbb{S}_i^{(2n-1,2n)}$ for R2-type parties, $U^\pi_i=\overrightarrow{\Pi}_i$ for T1-type parties, and $U^\pi_i=\overleftarrow{\Pi}_i$ for T2-type parties. In fact, the definitions of $U_i^\pi$ for R1-type and R2-type parties can be exchanged, and the definitions of $U_i^\pi$ for T1-type and T2-type parties can also be exchanged. This can be proved by the permutation invariant properties of $M_{2n}^\pi$.

Define $\Pi_{\mathrm{inv}}$ to be the $2n$-th order inversely permutation operator, which satisfies $\Pi_{\mathrm{inv}}[\rho^{(1)}\otimes\rho^{(2)}\cdots\otimes\rho^{(2n)}]\Pi_{\mathrm{inv}}^\dagger=\rho^{(2n)}\otimes\rho^{(2n-1)}\cdots\otimes\rho^{(1)}$. Because the $2n$ copies of $\rho$ in $M_{2n}^\pi$ are identical, we have
\begin{equation}
\begin{split}
M_{2n}^\pi=\tr\left[\left(\bigotimes_{i=1}^kU_i^\pi\right)\Pi_{\mathrm{inv}}\left(\rho^{\otimes 2n}\right)\Pi_{\mathrm{inv}}^\dagger\right]=\tr\left[\Pi_{\mathrm{inv}}^\dagger\left(\bigotimes_{i=1}^kU_i^\pi\right)\Pi_{\mathrm{inv}}\left(\rho^{\otimes 2n}\right)\right]=\tr\left[\left(\bigotimes_{i=1}^k\Pi_{\mathrm{inv},i}^\dagger U_i^\pi\Pi_{\mathrm{inv},i}\right)\left(\rho^{\otimes 2n}\right)\right].
\end{split}
\end{equation}
For R1-type and R2-type parties, $\Pi_{\mathrm{inv},i}^\dagger U_i^\pi\Pi_{\mathrm{inv},i}=U_i^\pi$, while for T1-type parties, $\Pi_{\mathrm{inv},i}^\dagger \overrightarrow{\Pi}_i\Pi_{\mathrm{inv},i}=\overrightarrow{\Pi}_i^\dagger=\overleftarrow{\Pi}_i$, and for T2-type parties, $\Pi_{\mathrm{inv},i}^\dagger \overleftarrow{\Pi}_i\Pi_{\mathrm{inv},i}=\overleftarrow{\Pi}_i^\dagger=\overrightarrow{\Pi}_i$. If we replace $\Pi_{\mathrm{inv},i}$ with $\overrightarrow{\Pi}$,
\begin{equation}
\begin{split}
M_{2n}^\pi=\tr\left[\left(\bigotimes_{i=1}^kU_i^\pi\right)\overrightarrow{\Pi}\left(\rho^{\otimes 2n}\right)\overrightarrow{\Pi}^\dagger\right]=\tr\left[\overrightarrow{\Pi}^\dagger\left(\bigotimes_{i=1}^kU_i^\pi\right)\overrightarrow{\Pi}\left(\rho^{\otimes 2n}\right)\right]=\tr\left[\left(\bigotimes_{i=1}^k\overrightarrow{\Pi}_i^\dagger U_i^\pi\overrightarrow{\Pi}_i\right)\left(\rho^{\otimes 2n}\right)\right].
\end{split}
\end{equation}
In this scenario, the definitions of $U_i^\pi$ of T1-type and T2-type parties will not change. While for R1-type parties, $\overrightarrow{\Pi}^\dagger_i[\mathbb{S}_i^{(2,3)}\otimes \cdots\otimes\mathbb{S}_i^{(2n,1)}]\overrightarrow{\Pi}_i=\mathbb{S}_i^{(1,2)}\otimes \cdots\otimes\mathbb{S}_i^{(2n-1,2n)}$, and for R2-type parties, $\overrightarrow{\Pi}^\dagger_i[\mathbb{S}_i^{(1,2)}\otimes \cdots\otimes\mathbb{S}_i^{(2n-1,2n)}]\overrightarrow{\Pi}_i=\mathbb{S}_i^{(2,3)}\otimes \cdots\otimes\mathbb{S}_i^{(2n,1)}$. Similarly, if we adopt $\overrightarrow{\Pi}\Pi_{\mathrm{inv}}$ to rotate these $2n$ copies of $\rho$, we can prove that the definitions of $U_i^\pi$ for R1-type and R2-type parties, and the definitions for T1-type and T2-type parties can be exchanged simultaneously.

\comments{
However, this is not in an observable form and is not measurable. We need to prove that this equation also holds when changing $U_i^\pi$ to $U_i^{\pi\dagger}$, which is equivalent to exchanging the definition of $U_i^\pi$ for T1-type and T2-type parties.

In fact, the definitions of $U_i^\pi$ for R1-type parties and R2-type parties can be exchanged, and so can T1-type parties and T2-type parties. This can be proved by changing the order of these $2n$ copies of $\rho$. Imagine in Eq.~\eqref{eq:M2nunitary}, all these $2n$ identical copies of $\rho$ have virtual label, $M_{2n}^\pi=\tr\left[\left(\bigotimes_{i=1}^kU_i^\pi\right)\rho^{(1)}\otimes\rho^{(2)}\otimes\cdots\otimes\rho^{(2n)}\right]$. If we inverse the order of these $2n$ copies of $\rho$ and ask the equality to hold for the new order, we need to change the definition of $U_i^\pi$ for different types of parties, $M_{2n}^\pi=\tr\left[\left(\bigotimes_{i=1}^kU_i^{\pi'}\right)\rho^{(2n)}\otimes\rho^{(2n-1)}\otimes\cdots\otimes\rho^{(1)}\right]$. Without loss of generality, we use the fourth order moment to show the change of $U_i^\pi$ for the four types of parties,
\begin{eqnarray}
\begin{tabular}{c}
\includegraphics[scale=0.14]{T1_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{T1_4_outorder.png}
\end{tabular}
\ , \
\begin{tabular}{c}
\includegraphics[scale=0.14]{T2_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{T2_4_outorder.png}
\end{tabular}
\ , \
\begin{tabular}{c}
\includegraphics[scale=0.14]{R1_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{R1_4_outorder.png}
\end{tabular}
\ , \
\begin{tabular}{c}
\includegraphics[scale=0.14]{R2_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{R2_4_outorder.png}
\end{tabular}
\end{eqnarray}
Therefore, if set $U_i^\pi=\overrightarrow{\Pi}^\dagger=\overleftarrow{\Pi}$ for T1-type parties, $U_i^\pi=\overleftarrow{\Pi}^\dagger=\overrightarrow{\Pi}$ for T2-type parties, and keep the definitions of $U_i^\pi$ for R1-type and R2-type parties, Eq.~\eqref{eq:M2nunitary} also holds. If we cyclically change the order of these $2n$ copies of $\rho$, $M_{2n}^\pi=\tr\left[\left(\bigotimes_{i=1}^kU_i^{\pi'}\right)\rho^{(2n)}\otimes\rho^{(1)}\otimes\cdots\otimes\rho^{(2n-1)}\right]$, the definitions of $U_i^\pi$ for T1-type parties and T2-type parties will not change, while the definitions of $U_i^\pi$ for R1-type parties and R2-type parties will exchange, which can be shown also by the fouth order moments,
\begin{eqnarray}
\begin{tabular}{c}
\includegraphics[scale=0.14]{T1_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{T1_4_perorder.png}
\end{tabular}
\ , \
\begin{tabular}{c}
\includegraphics[scale=0.14]{T2_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{T2_4_perorder.png}
\end{tabular}
\ , \
\begin{tabular}{c}
\includegraphics[scale=0.14]{R1_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{R1_4_perorder.png}
\end{tabular}
\ , \
\begin{tabular}{c}
\includegraphics[scale=0.14]{R2_4_inorder.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.14]{R2_4_perorder.png}
\end{tabular}
\end{eqnarray}
If we inverse and cyclically change the order of these $2n$ copies of $\rho$, $\rho^{(1)}\otimes\rho^{(2)}\otimes\cdots\otimes\rho^{(2n)}\to\rho^{(2n)}\otimes\rho^{(2n-1)}\otimes\cdots\otimes\rho^{(1)}\to\rho^{(1)}\otimes\rho^{(2n)}\otimes\cdots\otimes\rho^{(2)}$, we can prove that Eq.~\eqref{eq:M2nunitary} also holds when changing the definitions of $U_i^\pi$ for T1-type parties and T2-type parties, and the definitions of $U_i^\pi$ for R1-type parties and R2-type parties simultaneously.

Therefore, we can rewrite Eq.~\eqref{eq:M2nunitary} into observable form
\begin{equation}
\begin{split}
M_{2n}^\pi=\frac{1}{2}\tr\left[\left(\bigotimes_{i=1}^kU_i^\pi+h.c.\right)\rho^{\otimes 2n}\right],
\end{split}
\end{equation}
where $U_i^\pi$ follows the same definition in Eq.\eqref{eq:M2nunitary}.
\end{proof}
Take a four-partite state $ρ$ as an example. The four parties of $ρ$ are all different types with respect to $ℛ_π(·)$, T1-type, T2-type, R1-type, and R2-type respectively. It can be represented using tensor network as
\begin{eqnarray}
\mathcal{R}_\pi(\rho)=
\begin{tabular}{c}
\includegraphics[scale=0.12]{Rpirho.png}
\end{tabular}
\end{eqnarray}
Because $ℛ_π$ and $ℛ_π^†$ are mirror symmetric, the fourth moment, $M_4^π=[(ℛ_πℛ_π^†)^2]$ can be graphically represented as
\begin{eqnarray}
\begin{tabular}{c}
\includegraphics[scale=0.12]{M4intrace.png}
\end{tabular}
\right)
\begin{tabular}{c}
\includegraphics[scale=0.12]{M4outtrace.png}
\end{tabular}
\end{eqnarray}
As we said, one can find that the connection of legs belonging to T-type parties is among the four copies of $ρ$, while the connection of R-type legs is between two neighboring $ρ$. Hence, $M_4^π$ is equivalent to measuring
\begin{equation}
\begin{split}
\end{split}
\end{equation}
on $4$ copies of $ρ$.

\subsection{Optimization Problem}\label{subsec:optimization}
\begin{theorem}\label{theorem:Lagrange Multipler}
The minimum value of $\norm{\mathcal{R}_\pi}$ given $M_2^\pi$, ..., $M_{2n}^\pi$ is reached when there are at most $n$ non-zero $\lambda_i$s. Thus, denote the solution of this problem to be $E_{2n}^\pi(\rho)$, the minimum value of $\sum_i\lambda_i$ is the solution of the following optimization problem,
\begin{equation}
\begin{split}
\min_{q_1,\cdots,q_n\in\mathbb{N}}& E_{2n}^\pi(\rho)= q_1\lambda_1+q_2\lambda_2+\cdots+q_n\lambda_n\\
\mathrm{s.t.}& \ \sum_{i=1}^nq_i\lambda_i^2=M_2^\pi,\cdots, \ \sum_{i=1}^nq_i\lambda_i^{2n}=M_{2n}^\pi\\
& \ q_1+q_2+\cdots+q_n\le L
\end{split}
\end{equation}
\end{theorem}
\begin{proof}
The original optimization problem:
\begin{equation}\label{eq:biopt}
\begin{split}
\min_{\{\lambda_i\}} \ \ &\sum_{i=1}^L \lambda_i\\
\text{subject to} \ \ &\sum_{i=1}^L \lambda_i^{2}=M_{2}^\pi \\
&\sum_{i=1}^L \lambda_i^{2n}=M_{2n}^\pi \\
& \lambda_1\ge \lambda_2\ge\cdots\ge \lambda_L\ge 0.
\end{split}
\end{equation}
Here, $L$ is the number of the singular values. Let $\lambda_i=x_L^2+\cdots+x_i^2$ with $x_1,\cdots,x_L\in \mathbb{R}.$ The above optimization problem can be written as
\begin{equation}
\begin{split}
\min_{\{x_i\}} \ \ &\sum_{i=1}^L\left( x_L^2+\cdots+x_i^2\right)\\
\text{subject to} \ \ &\sum_{i=1}^L (x_L^2+\cdots+x_i^2)^{2}=M_{2}^\pi \\
&\sum_{i=1}^L (x_L^2+\cdots+x_i^2)^{2n}=M_{2n}^\pi \\
&x_1,\cdots,x_L\in \mathbb{R}.
\end{split}
\end{equation}
Define the Lagrange function as
\begin{equation}
\begin{split}
\mathcal{L}&=\sum_{i=1}^L  (x_L^2+\cdots+x_i^2)+\alpha_1\Big(\sum_{i=1}^L (x_L^2+\cdots+x_i^2)^{2}-M_{2}^\pi\Big)\\
&+\cdots +\alpha_n\left(\sum_{i=1}^L (x_L^2+\cdots+x_i^2)^{2n}-M_{2n}^\pi\right).
\end{split}
\end{equation}
Take partial derivative with respect to all variables and the minimum value in the optimization problem can be achieved at the points where all derivatives are equal to $0.$ Therefore, we have
\begin{equation}
\begin{split}
0=&\frac{\partial \mathcal{L}}{\partial x_L } =2L x_L+4\alpha_1 x_L\sum_{i=1}^L (x_L^2+\cdots+x_i^2)+\cdots+4n\alpha_n x_L\sum_{i=1}^L (x_L^2+\cdots+x_i^2)^{2n-1},\\
0=&\frac{\partial \mathcal{L}}{\partial x_{1} }=2x_1+4\alpha_1 x_1 (x_L^2+\cdots+x_1^2)+\cdots+4n\alpha_n x_1(x_L^2+\cdots+x_1^2)^{2n-1},
\end{split}
\end{equation}
substituting $\lambda_i=x_L^2+\cdots+x_i^2$ into this equation,
\begin{equation}
\begin{split}
\end{split}
\end{equation}
These above equations indicate that the extreme points satisfy
\begin{equation}\label{eq:origin solution}
\begin{split}
\lambda_L=&0\quad \text{or}  \quad L+2\alpha_1\sum_{i=1}^L \lambda_i+\cdots+2n\alpha_n\sum_{i=1}^L \lambda_i^{2n-1}=0,\\
\lambda_L=&\lambda_{L-1}\quad \text{or}  \quad L-1+2\alpha_1\sum_{i=1}^{L-1} \lambda_i+\cdots+2n\alpha_n\sum_{i=1}^{L-1} \lambda_i^{2n-1}=0,\\
\lambda_2=&\lambda_{1}\quad \text{or}  \quad 1+2\alpha_1\lambda_1+\cdots+2n\alpha_n\lambda_1^{2n-1}=0.
\end{split}
\end{equation}
After simple analysis, there's only two possible kinds of $\lambda_i$ that satisfies the above equations: $\lambda_i=0$ or are the roots of the high-degree equation
\begin{equation}\label{eq:extreme}
\end{equation}
Notice that all even-degree terms are zero and
\begin{equation}
\end{equation}
is an odd function. As a consequence, Eq.~\eqref{eq:extreme} has at most $n$ positive roots. This implies that the set $\{\lambda_i\}$ have at most $n$ different positive values and the rest are $0$. Assume the $n$ different positive values are $\lambda_1\ge \cdots\ge \lambda_n\ge 0$. Denote the number of $\lambda_i$ as $q_i$ and the sum satisfies $q_1+\cdots +q_n\le L$. The optimization problem can be reduced to a new simpler optimization problem.

\begin{equation}
\begin{split}
\min_{\{q_i\in\mathbb{N}\}} \ \ &q_1 \lambda_1+\cdots+q_n \lambda_n\\
\text{subject to} \ \ &\sum_{i=1}^n q_i\lambda_i^{2}=M_{2}^\pi \\
&\sum_{i=1}^n q_i\lambda_i^{2n}=M_{2n}^\pi \\
& \lambda_1\ge \lambda_2\ge\cdots\ge \lambda_n\ge 0\\
& q_1+\cdots q_n\le L.
\end{split}
\end{equation}
\end{proof}

Particularly, we give an analytical solution in the case $n=2$. The optimization problem is
\begin{equation}\label{eq:optimizationM2M4}
\begin{split}
\min_{\{\lambda_i\}} \ \ &\sum_{i=1}^L \lambda_i\\
\text{subject to} \ \ &\sum_{i=1}^L \lambda_i^{2}=M_{2}^\pi \\
&\sum_{i=1}^L \lambda_i^{4}=M_{4}^\pi \\
& \lambda_1\ge \lambda_2\ge\cdots\ge \lambda_L\ge 0.
\end{split}
\end{equation}
To solve this problem, we need to adopt a theorem in [\citenamefont{Berry and Sanders}, 2003].
\begin{fact}
Suppose $\{p_1,p_2,\cdots,p_{L}\}$ is a probability distribution, namely,  $p_i\geq 0$ and $\sum_i p_i=1.$ The sum $H_f=\sum_i f(p_i)$ is called an entropy measure, if $f(\cdot)$ is a function satisfying the following conditions,
\begin{enumerate}
\item f(0)=0,
\item $f(\cdot)$ is strictly convex or strictly concave,
\item the first derivative $f^{\prime}$ exists and is continuous in the interval $(0,1)$.
\end{enumerate}

Let $f(\cdot)$ and $g(\cdot)$ be two functions satisfying the above conditions. Thus, $H_f=\sum_i f(p_i)$ and $H_g=\sum_i g(p_i)$ are two entropy measures. Define $\tilde{f}^{\prime}(g^{\prime})=f^{\prime}(p(g^{\prime}))$, which is formulated by treating $f^{\prime}=\frac{df}{dp}$ as the function of $g^{\prime}=\frac{dg}{dp}$. If $\tilde{f}^{\prime}(g^{\prime})$ is a strictly convex function of $g^{\prime}$, then the solution of the optimization problem
\begin{equation}\label{eq:optimizationentropy}
\begin{split}
\min_{\{p_i\ge 0\}} \ \ &H_f=\sum_{i=1}^L f(p_i)\\
subject \ to \ \ &\sum_{i=1}^L p_i=1\\
& \sum_{i=1}^Lg(p_i)=H_g
\end{split}
\end{equation}
is obtained when
\begin{equation}
\begin{split}
\{p_i\}=\{p_1,\cdots,p_1,p_2,0,\cdots,0\},\qquad p_2=1-p_1\left \lfloor \frac{1}{p_1}  \right \rfloor< p_1.
\end{split}
\end{equation}
If $\tilde{f}^{\prime}(g^{\prime})$ is a strictly concave function of $g^{\prime}$, the minimum is obtained when
\begin{equation}
\{p_i\}=\{p_1,p_2,\cdots,p_2\},\qquad p_2=\frac{1-p_1}{d-1}\leq p_1.
\end{equation}

\end{fact}
Substituting $p_i=λ_i^2$ , $f(p)=√(p)$  and $g(p)=p^2$ and ignore the normalization constant, Eq.~\eqref{eq:optimizationentropy} can be converted to our target optimization problem, Eq.~\eqref{eq:optimizationM2M4}. It is easy to verify that the two functions, $f(p)=√(p)$  and $g(p)=p^2$, satisfy the required three conditions, and $f̃^'(g^')=1/(√(2g^'))$ is strictly convex. Therefore, the minimum value is achieved when
\begin{equation}
\{\lambda_i\}=\{\lambda_1,\cdots,\lambda_1,\lambda_2,0,\cdots,0\}.
\end{equation}
Suppose the number of $λ_1$ is $q$, then institute it into the constraint conditions, one gets
\begin{equation}
\begin{split}
\lambda_1 &=\frac{1}{q\sqrt{q+1}}\sqrt{q(qM_2^\pi+U)}
\lambda_2 &=\frac{1}{\sqrt{q+1}}\sqrt{M_2^\pi-U},
\end{split}
\end{equation}
with $U=√(q(q+1)M_4^π-q(M_2^π)^2)$. To make this a valid solution, there is only one possible value of $q$, which is $q=⌊(M_2^π)^2/M_4 ⌋$. Therefore, the solution of the target optimization problem Eq.~\eqref{eq:optimizationM2M4} is
\begin{equation}
\begin{split}
\end{split}
\end{equation}

\subsection{Entanglement Structure Detection}\label{sec:structure}
In the multipartite system, when relaxing the restriction of LOCC, states can have much more complicated entanglement structures [\citenamefont{Lu et~al.}, 2018] than in the bipartite scenario. Concepts like entanglement depth and intactness naturally arise. A $k$-partite state $ρ$ is called $t$-separable iff it can be written as
\begin{equation}
\begin{split}
\rho=\sum_{g_1\cup g_2\cup\cdots\cup g_t=[k]}\sum_i p_{\{g\},i}\rho_{g_1}^i\otimes \rho_{g_2}^i\otimes\cdots\otimes\rho_{g_t}^i,
\end{split}
\end{equation}
where $g_1$, ..., $g_t$ are $t$ disjoint non-empty sets, ${p_{g},i}$ are probabilities satisfying the normalization condition, $∑_{g},ip_{g},i=1$. The largest number of $t$ for a given $k$-partite state $ρ$ is called the \emph{entanglement intactness} of $ρ$.

To generalize our moment-based permutation criteria to entanglement structure detection, we need to first develop the norm-based permutation entanglement structure criteria, and then prove that they can also be estimated by the permutation moments. Adopting the Ky Fan matrix norm [\citenamefont{Horn et~al.}, 1994], we find two indicators for entanglement intactness.

\begin{theorem}
Given a $k$-partite state $\rho$ with intactness $t$, it satisfies
\begin{equation}
\begin{split}
G_R(\rho)&=\sum_{g\subsetneq [k]}\norm{\mathcal{R}_{(2,3)}(\rho_{g,\bar{g}})}_{d^2}\le(2^k-2^t)d+(2^t-2),\\
G_T(\rho)&=\sum_{g\subsetneq [k]}\norm{\mathcal{R}_{(1,2)}(\rho_{g,\bar{g}})}_{d^2}\le(2^k-2^t)d+(2^t-2),
\end{split}
\end{equation}
where $\norm{\cdot}_{d^2}$ is the Ky Fan $d^2$ matrix norm, defined by the sum of $d^2$ largest singular values of the matrix; $d$ is an integer that is not greater than the dimension of the smallest party; $\rho_{g,\bar{g}}$ is the bipartite state constructed by treating $g$ and $\bar{g}$ as the two parties of $\rho$, where $\bar{g}$ is the complement of $g$; $\mathcal{R}_{(2,3)}$ and $\mathcal{R}_{(1,2)}$ are the bipartite index permutation operations acting on the indices of parties $g$ and $\bar{g}$.
\end{theorem}

When $k=3$, this criterion gives the genuine tripartite entanglement criterion proposed in Ref.~[\citenamefont{Li et~al.}, 2017].

\begin{proof}
Here, we use $G_R(\rho)$ as an example to prove our theorem, the proof for $G_T(\rho)$ is quite similar. According to the convexity of Ky Fan norm, these two functions are all convex:
\begin{equation}
\begin{split}
G_R(a\rho_1+b\rho_2)&=\sum_{g\subsetneq [k]}\norm{\mathcal{R}_{(2,3)}\left[a(\rho_1)_{g,\bar{g}}+b(\rho_2)_{g,\bar{g}}\right]}_{d^2}\\
&\le\sum_{g\subsetneq [k]}a\norm{\mathcal{R}_{(2,3)}\left[(\rho_1)_{g,\bar{g}}\right]}_{d^2}+b\norm{\mathcal{R}_{(2,3)}\left[(\rho_2)_{g,\bar{g}}\right]}_{d^2}\\
&=a G_R(\rho_1)+b G_R(\rho_2).
\end{split}
\end{equation}
Hence, for a state with intactness $t$, $\rho=\sum_{i,\psi_i\in t-int}p_i\ketbra{\psi_i}{\psi_i}$,
\begin{equation}
\begin{split}
G_R(\rho)\le \max_{\psi\in t-int}G_R(\psi),
\end{split}
\end{equation}
where we use $t-int$ to denote the set of states that have intactness $t$. For pure state with intactness $t$, $\ket{\psi}=\bigotimes_{i=1}^t\ket{\psi_i}$, there's $2^t-2$ non-trivial partition $g|\bar{g}$ that ensure each $\ket{\psi_i}$ is in either $g$ or $\bar{g}$. In this case, $g$ and $\bar{g}$ are two separate parts, $\norm{\mathcal{R}_{(2,3)}(\psi_{g,\bar{g}})}_{d^2}=1$ because $\tr\left\{[\mathcal{R}_{(2,3)}(\psi_{g,\bar{g}})\mathcal{R}_{(2,3)}(\psi_{g,\bar{g}})^\dagger]^n\right\}=1$ for all $n$. In other cases, there exists some $\psi_i$ that distributed both in $g$ and $\bar{g}$, hence $g$ and $\bar{g}$ are two entangled parties. Also use singular values to denote the norm, $\norm{\mathcal{R}_{(2,3)}(\psi_{g,\bar{g}})}_{d^2}=\sum_{i=1}^{d^2}\lambda_i$, where $\lambda_i$ are arranged in decreasing order, the purity condition $\tr(\psi^2)=1$ gives a constrain that $\sum_{i=1}^{d^2}\lambda_i^2\le 1$. Recalling that singular values are all non-negative, hence the maximum value of $\sum_{i=1}^{d^2}\lambda_i$ under this purity constrain is $d$.

Considering other constraints in addition to purity constraint, the maximum will be less than $d$. In fact, the exact maximum value depends on the number of $\psi_i$ that is not fully contained in $g$ or $\bar{g}$. We will leave the calculation of the exact upper bound in our future work. Therefore, $(2^k-2^t)d+(2^t-2)$ generally cannot be achieved by some $k$-partite states with intactness $t$. However, this bound is nontrivial because $k$-partite GHZ state reaches $(2^k-2)d$ which is larger than $(2^k-2^t)d+(2^t-2)$ for any $t\ge 2$.

\end{proof}

The lower bound of Ky Fan norm can also be constructed using moments. Here we only consider the case knowing the second and fourth moments. The optimization problem for entanglement structure detection is similar as one in Sec.~\ref{subsec:optimization}, and can be simplified as
\begin{equation}\label{eq:gmeopt}
\begin{split}
\min_{\lambda_i} \ \ &\sum_{i=1}^l \lambda_i\\
\text{subject to} \ \ &\sum_{i=1}^L \lambda_i^{2}=M_{2} \\
&\sum_{i=1}^L \lambda_i^{4}=M_{4} \\
& \lambda_1\ge \lambda_2\ge\cdots\ge \lambda_L\ge 0\\
& l\le L.
\end{split}
\end{equation}
Using the similar strategy to simplify this optimization problem, it can be obtained that the extreme points satisfy
\begin{equation}
\begin{split}
\lambda_L&=0\quad \text{or}  \quad l+2\alpha\sum_{i=1}^L \lambda_i+4\beta\sum_{i=1}^L \lambda_i^{3}=0,\\
\lambda_{l+1}&=\lambda_{l}\quad \text{or}  \quad l+2\alpha\sum_{i=1}^l \lambda_i+4\beta\sum_{i=1}^l \lambda_i^{3}=0,\\
\lambda_l&=\lambda_{l-1}\quad \text{or}  \quad l-1+2\alpha\sum_{i=1}^{l-1} \lambda_i+4\beta\sum_{i=1}^{l-1} \lambda_i^{3}=0,\\
\lambda_2&=\lambda_{1}\quad \text{or}  \quad 1+2\alpha \lambda_1+4\beta\lambda_1^{3}=0,
\end{split}
\end{equation}
where $α$ and $β$ are Lagrange coefficients. Compared with the original optimization problem, Eq.~\eqref{eq:origin solution}, the above equations are different in constant terms of the first $(L-l)$ lines. We will see in the following analysis that this difference will make the problem much harder than the original problem. We discuss the solution in different cases.

\begin{figure}
\centering
\includegraphics[width=15cm]{gmeoptimization.pdf}
\caption{Different cases in entanglement structure optimization problem.}
\label{fig:gmecase}
\end{figure}

1).  $λ_l≠λ_l+1$.

For $i≤l,$ $λ_i$ is the root of the equation $4βλ^3+2αλ=-1$.  For $i>l$,  $λ_i$ is the root of the equation $4βλ^3+2αλ=0$ and obviously, $0$ is one root. For the convenience of discussion, we draw the distribution diagrams of roots Fig.~\ref{fig:gmecase}a1) and Fig.~\ref{fig:gmecase}a2), and there are two possible cases, depending on the sign of $α$ and $β$. For the first case, as shown in  Fig~.\ref{fig:gmecase}a1), $λ_i$ must equal to 0 for all $i>l$ because the possible nonzero solution is greater than $λ_1$ and $λ_2$. Hence, the original optimization problem becomes a much simpler minimization problem,
\begin{equation}
\begin{split}
\min \ \ &f=q_1\lambda_1+q_2\lambda_2 \\
\text{subject to} \ \ &q_1\lambda_1^2+q_2\lambda_2^2=M_2\\
&q_1\lambda_1^4+q_2\lambda_2^4 =M_4,\\
&\lambda_1\ge \lambda_2> 0,\quad q_1,q_2\in\mathbb{N}, \quad q_1+q_2= l.
\end{split}
\end{equation}
For the second case, as shown in Fig~.\ref{fig:gmecase}a2), $λ_i$ only has one possible value $λ_1$ for all $i≤l$, while has a possible positive value $λ_2$ for $i>l$. Hence, the original optimization problem becomes,
\begin{equation}
\begin{split}
\min \ \ &f=l\lambda_1+q_2\lambda_2 \\
\text{subject to} \ \ &l\lambda_1^2+q_2\lambda_2^2=M_2\\
&l\lambda_1^4+q_2\lambda_2^4 =M_4,\\
&\lambda_1\ge \lambda_2> 0,\quad q_2\in\mathbb{N},\\
&q_2\leq L-l.
\end{split}
\end{equation}
In fact, it is possible to find the analytical solutions for these two cases, by the ordinary derivation process. However, in the following cases, finding analytical solutions is not likely to be done.

2).  $λ_l= λ_l+1$ and there exists an positive integer $u$ such that $u<l$, $λ_u≠λ_u+1$, $λ_u+1=λ_u+2=⋯=λ_L=0$.
Then $λ_1,⋯,λ_u$ are roots of the equation $4βλ^3+2αλ+1 =0$, the distribution is shown in Fig.~\ref{fig:gmecase}a2) which is quite similar to the one shown in Fig.~\ref{fig:gmecase}a1), and can be unified with that case to get a more general minimization problem:
\begin{equation}\label{eq:gmecase1}
\begin{split}
\min \ \ &f=q_1\lambda_1+q_2\lambda_2 \\
\text{subject to} \ \ &q_1\lambda_1^2+q_2\lambda_2^2=M_2\\
&q_1\lambda_1^4+q_2\lambda_2^4 =M_4,\\
&\lambda_1\ge \lambda_2> 0,\quad q_1,q_2\in\mathbb{N},\\
&q_1+q_2\le l.
\end{split}
\end{equation}

3). If  $λ_1=⋯=λ_l+1=⋯=λ_q≠λ_q+1$ and $q≥l+1$, then $λ_1,⋯,λ_q$ are roots of the equation $4βλ^3+2αλ+l/q=0$ and $λ_q+1,⋯,λ_L$ are roots of the equation $4βλ^3+2αλ=0$. As shown in Fig~.\ref{fig:gmecase}b2), there is only one possible case.  This case can be unified with the case shown in Fig.~\ref{fig:gmecase}b1), to give a more general optimization problem:
\begin{equation}\label{eq:gmecase2}
\begin{split}
\min \ \ &f=q_1\lambda_1+q_2\lambda_2 \\
\text{subject to} \ \ &q_1\lambda_1^2+q_2\lambda_2^2=M_2\\
&q_1\lambda_1^4+q_2\lambda_2^4 =M_4,\\
&\lambda_1\ge \lambda_2> 0,\quad,q_1,q_2\in\mathbb{N},\\
&l\le q_1,\quad q_1+q_2\le L.
\end{split}
\end{equation}

4).  $λ_l= λ_l+1$ and there exists integers $u$ and $v$ satisfying $0<u<l$ and $v>l+1$, such that $λ_u≠λ_u+1=⋯=λ_l=λ_l+1=⋯=λ_v≠λ_v+1=0$. Then $λ_1,⋯,λ_u$ are roots of the equation $4βλ^3+2αλ=-1$ and $λ_u+1,⋯,λ_v$ are roots of the equation $4βλ^3+2αλ=-l-u/v-u$. The last terms $λ_v+1,⋯,λ_L$ are roots of $4βλ^3+2αλ=0$. As shown in Fig~.\ref{fig:gmecase}c) and Fig.~\ref{fig:gmecase}d), there are also two cases dependent on the sign of $α$ and $β$. Now we discuss them in details.

For the first case, as shown in Fig~.\ref{fig:gmecase}c), $λ_i$ has two possible values for $i≤u$ and one possible value for $u<i≤v$. In this scenario, the optimization problem becomes
\begin{equation}
\begin{split}
\min \ \ &f=q_1\lambda_1+q_2\lambda_2 +q_3\lambda_3\\
\text{subject to} \ \ &q_1\lambda_1^2+q_2\lambda_2^2+q_3\lambda_3^2=M_2\\
&q_1\lambda_1^4+q_2\lambda_2^4 +q_3\lambda_3^4=M_4,\\
&4\beta \lambda_1^3+2\alpha \lambda_1+1=0,\\
&4\beta \lambda_2^3+2\alpha \lambda_2+1=0,\\
&4\beta \lambda_3^3+2\alpha \lambda_3+\frac{l-q_1-q_2}{q_3}=0,\\
&\lambda_1\ge \lambda_2\geq\lambda_3 >0,\quad q_1,q_2,q_3\in\mathbb{N},\\
&\quad q_1+q_2<l,\quad l<q_1+q_2+q_3\le L.
\end{split}
\end{equation}
This optimization problem can be further simplified to
\begin{equation}\label{eq:gmecase3}
\begin{split}
\min \ \ &f=q_1\lambda_1+q_2\lambda_2 +q_3\lambda_3\\
\text{subject to} \ \ &q_1\lambda_1^2+q_2\lambda_2^2+q_3\lambda_3^2=M_2\\
&q_1\lambda_1^4+q_2\lambda_2^4 +q_3\lambda_3^4=M_4,\\
&\lambda_1\ge \lambda_2\geq\lambda_3 >0,\quad q_1,q_2,q_3\in\mathbb{N},\\
&\quad q_1+q_2<l,\quad l<q_1+q_2+q_3\le L.
\end{split}
\end{equation}
For the second case, as seen in Fig~.\ref{fig:gmecase}d), $λ_i$ has one possible value in each of the three intervals, $1≤i≤u$, $u+1≤i≤v$ and $v+1≤i≤L$. Then the optimization problem becomes
\begin{equation}\label{eq:gmecase4}
\begin{split}
\min \ \ &f=q_1\lambda_1+q_2\lambda_2 +q_3\lambda_3\\
\text{subject to} \ \ &q_1\lambda_1^2+q_2\lambda_2^2+q_3\lambda_3^2=M_2\\
&q_1\lambda_1^4+q_2\lambda_2^4 +q_3\lambda_3^4=M_4,\\
&\lambda_1\ge \lambda_2\geq\lambda_3 >0,\quad q_1,q_2,q_3\in\mathbb{N},\\
& q_1<l , l< q_1+q_2<L, q_1+q_2+q_3\le L.
\end{split}
\end{equation}

Therefore, to solve the entanglement structure optimization problem, Eq.~\eqref{eq:gmeopt}, what we need to do is to find the minimum result of the four much easier optimization problems, Eq.~\eqref{eq:gmecase1}, Eq.~\eqref{eq:gmecase2}, Eq.~\eqref{eq:gmecase3} and Eq.~\eqref{eq:gmecase4}. Although it is also hard to find the analytical result, the complexity of this problem has been greatly reduced.

\section{Measurement Protocols of Permutation Moments and The Statistical Analysis}\label{sec:protocols}
\subsection{Shadow Estimation and Statistical Analysis}\label{subsec:shadow}
Shadow tomography [\citenamefont{Huang et~al.}, 2020] is a systematic framework which can help us to benchmark the property of an unknown quantum system with only partial knowledge. Compared with former protocols, shadow tomography is much more resource-saving than full tomography and it does not need joint operations among multiple copies to estimate nonlinear functions.  Shadow tomography consists of two phases: quantum measurement and data postprocessing. In quantum measurement phase, we adopt $U$ to rotate $N$-qubit state $ρ$ to $UρU^†$, and measure it in computational basis ${|s⃗⟩ =|s_1,⋯,s_N⟩}$. The measurement results and the rotation unitary together construct unbiased estimators of $ρ$:
\begin{equation}\label{eq:global shadow}
\begin{split}
\hat{\rho}=(2^N+1)U^\dagger\ketbra{\vec{s}}{\vec{s}}U-\mathbb{I}_{2^N}
\end{split}
\end{equation}
for $U$ randomly chosen from $N$-qubit Clifford group, and
\begin{equation}\label{eq:local shadow}
\begin{split}
\hat{\rho}=\bigotimes_{i=1}^N\left(3u_i^\dagger\ketbra{s_i}{s_i}u_i-\mathbb{I}_2\right)
\end{split}
\end{equation}
for $U=⊗_i=1^Nu_i$ with each $u_i$ randomly chosen from single qubit Clifford group.

Recalling that $ℛ_π$ is just the index rearrangement of $ρ$, so shadow tomography can also construct an unbiased estimator of it
\begin{equation}
\begin{split}
\hat{\mathcal{R}}_\pi=\mathcal{R}_\pi(\hat{\rho}), \mathbb{E}[\hat{\mathcal{R}}_\pi]=\mathcal{R}_\pi(\mathbb{E}\hat{\rho})=\mathcal{R}_\pi.
\end{split}
\end{equation}
Hence, the unbiased estimator for $M_2n^π$ can be constructed using these $ℛ̂_π$
\begin{equation}
\begin{split}
\hat{M}_{2n}^\pi=\frac{(M-2n)!}{M!}\sum_{\{i_1,\cdots, i_{2n}\}\subset[M]}\sum_{\sigma\in\mathcal{S}_{2n}}\tr[\hat{\mathcal{R}}_{\pi\sigma(i_1)}\cdots\hat{\mathcal{R}}_{\pi\sigma(i_{2n})}^\dagger],
\end{split}
\end{equation}
where $M$ is the number of the $ℛ̂_π$ we prepared using shadow tomography, $𝒮_2n$ is the 2n-order permutation group, and $ℛ̂_πσ(i)$ is the $σ(i)$-th copy of $ℛ̂_π$. The unbiaseness of $M̂_2n$ is easy to be proved:
\begin{equation}
\begin{split}
\mbb{E}\left(\hat{M}_{2n}^\pi\right)&=\frac{(M-2n)!}{M!}\sum_{\{i_1,\cdots, i_{2n}\}\subset[M]}\sum_{\sigma\in\mathcal{S}_{2n}}\tr[\mbb{E}\left(\hat{\mathcal{R}}_{\pi\sigma(i_1)}\right)\cdots\mbb{E}\left(\hat{\mathcal{R}}_{\pi\sigma(i_{2n})}^\dagger\right)]\\
&=\frac{(M-2n)!}{M!}\sum_{\{i_1,\cdots, i_{2n}\}\subset[M]}\sum_{\sigma\in\mathcal{S}_{2n}}\tr\left[\mathcal{R}_\pi\cdots\mathcal{R}_\pi^\dagger\right]\\
\end{split}
\end{equation}

In fact, there exists many other protocols to measure nonlinear functions of state, like randomized measurements[\citenamefont{Brydges et~al.}, 2019, \citenamefont{Elben et~al.}, 2019], which has been proved to be more efficient than shadow estimation in some cases [\citenamefont{Elben et~al.}, 2020, \citenamefont{Zhou et~al.}, 2020, \citenamefont{Liu et~al.}, 2022]. However, unlike shadow estimation, randomized measurements cannot always be easily conducted by qubit level operation, especially for observables like higher order permutation operators. Therefore, except for some special cases we will discuss in Sec.~\ref{subsec:random}, shadow estimation is a practical choice to measure the permutation moments.

In the following context, we take $M_4^(2,3)(ρ_AB)$ as an example to illustrate the variance of different protocols, which is closely related to the bipartite entanglement criterion proposed in this work, and the conclusion can be easily generalized to other permutation moments. To analyze the error scaling of estimators constructed by shadow estimation, here we adopt the conclusions in [\citenamefont{Huang et~al.}, 2020]:
\begin{fact}\label{fact:shadow variance}
The variance of predicting $\hat{O}=\tr(O\hat{\rho})$ on locally constructed $\hat{\rho}$, Eq.~\eqref{eq:local shadow}, scales like
\begin{equation}
\begin{split}
\mathrm{Var}(\hat{O})\leq 2^{\mathrm{locality}(\hat{O})}\tr(\hat{O}^2).
\end{split}
\end{equation}
For globally constructed $\hat{\rho}$, Eq.~\eqref{eq:global shadow}, the variance scales like
\begin{equation}
\begin{split}
\mathrm{Var}(\hat{O})\leq 3\tr(\hat{O}^2).
\end{split}
\end{equation}
\end{fact}

Given $M$ measurement results of shadow estimation, the unbiased estimator we constructed to predict $M_4^(2,3)$ is
\begin{equation}\label{eq:M4def}
\begin{split}
\hat{M}_4^{(2,3)}=\frac{4(M-4)!}{M!}\sum_{\{i,j,k,l\}\subset[M]}\sum_{\sigma\in\mathcal{S}_3}\tr\left[\hat{R}_{i}\hat{R}_{\sigma(j)}^\dagger\hat{R}_{\sigma(k)}\hat{R}_{\sigma(l)}^\dagger\right],
\end{split}
\end{equation}
where for simplicity, we use $R̂$ to represent $ℛ̂_(2,3)(ρ_AB)$. Here, because of the properties of trace function that it is invariant under permutation and taking Hermite conjugate, $[R̂_jR̂_i^†R̂_kR̂_l^†]=[R̂_kR̂_l^†R̂_jR̂_i^†]=[R̂_iR̂_j^†R̂_lR̂_k^†]$, we can always set $R̂_i$ in the first place. Combining Fact.~\ref{fact:shadow variance} and the definition Eq.~\eqref{eq:M4def}, one can prove that:

\begin{proposition}\label{proposition:M4shadowvariance}
According to Chebyshev's equation, to make sure our estimation of $M_4$ is accurately enough, which is to say, the estimator defined in Eq.\eqref{eq:M4def} satisfies $|\hat{M}_4^{(2,3)}-M_4^{(2,3)}|\le\epsilon$ with probability at least $1-\delta$, the number of snapshots needed is
\begin{equation}
\begin{split}
\end{split}
\end{equation}
for global shadow protocol, and
\begin{equation}
\begin{split}
\end{split}
\end{equation}
for local shadow protocol, where $D$ stands for the system dimension.
\end{proposition}

\begin{proof}
Following the definition of variance, we have
\begin{equation}
\begin{split}
\mathrm{Var}(\hat{M}_4^{(2,3)})=\mathbb{E}\left[\left(\hat{M}_4^{(2,3)}\right)^2\right]-\mathbb{E}\left[\hat{M}_4^{(2,3)}\right]^2=\mathbb{E}\left[\left(\hat{M}_4^{(2,3)}\right)^2\right]-\left(M_4^{(2,3)}\right)^2.
\end{split}
\end{equation}
substituting Eq.~\eqref{eq:M4def}, we get
\begin{equation}
\begin{split}
\left(\hat{M}_4^{(2,3)}\right)^2=\left(\frac{4(M-4)!}{M!}\right)^2\sum_{\{i,j,k,l\}\subset[M]}\sum_{\{i',j',k',l'\}\subset[M]}\left\{\sum_{\sigma\in\mathcal{S}_3}\tr\left[\hat{R}_{i}\hat{R}_{\sigma(j)}^\dagger\hat{R}_{\sigma(k)}\hat{R}_{\sigma(l)}^\dagger\right]\right\}\left\{\sum_{\pi\in\mathcal{S}_3}\tr\left[\hat{R}_{i'}\hat{R}_{\pi(j')}^\dagger\hat{R}_{\pi(k')}\hat{R}_{\pi(l')}^\dagger\right]\right\}
\end{split}
\end{equation}
To benefit our calculation, we can divide the summation terms in the R.H.S of the last equation into five groups, according to the number of equal indices between $\{i,j,k,l\}$ and $\{i',j',k',l'\}$. If there are $0\leq c\leq 4$ pairs of same indices, the total number of such terms are $\binom{M}{4}\binom{4}{c}\binom{M-4}{4-c}$, hence
\begin{equation}
\begin{split}
\mathbb{E}\left[\left(\hat{M}_4^{(2,3)}\right)^2\right]&=\left(\frac{4(M-4)!}{M!}\right)^2\binom{M}{4}\binom{4}{0}\binom{M-4}{4}36\tr\left(RR^\dagger RR^\dagger\right)^2\\
&+\left(\frac{4(M-4)!}{M!}\right)^2\binom{M}{4}\binom{4}{1}\binom{M-4}{3}36\mathbb{E}\left[\tr\left(RR^\dagger R\hat{R}^\dagger\right)^2\right]\\
\end{split}
\end{equation}
where $R=\mathcal{R}_{(2,3)}(\rho_{AB})$, and $\hat{R}_1$, $\hat{R}_2$, $\hat{R}_3$, and $\hat{R}_4$ are independent unbiased snapshots of $R$. Hence,
\begin{equation}\label{eq:M_4variance}
\begin{split}
\mathrm{Var}\left[\left(\hat{M}_4^{(2,3)}\right)\right]&=36\left(\frac{4(M-4)!}{M!}\right)^2\binom{M}{4}\binom{4}{1}\binom{M-4}{3}\mathrm{Var}\left[\tr(RR^\dagger R\hat{R}^\dagger)\right]\\
&\leq \frac{C_1}{M}\mathrm{Var}\left[\tr(RR^\dagger R\hat{R}^\dagger)\right]+\frac{C_2}{M^2}\mathrm{Var}\left[\tr(RR^\dagger\hat{R}_1\hat{R}^\dagger_2)+\tr(RR^\dagger\hat{R}_2\hat{R}^\dagger_1)+\tr(R\hat{R}^\dagger_1R\hat{R}^\dagger_2)\right]\\
\end{split}
\end{equation}
Here, $C_1$, $C_2$, $C_3$ and $C_4$ are constants independent of $M$. To apply Fact \ref{fact:shadow variance}, we need to figure out the exact form of $\hat{O}_1$, $\hat{O}_2$, $\hat{O}_3$ and $\hat{O}_4$. Denote $\hat{O}$ to be the target observable $\hat{O}=\mathbb{S}_A^{(1,2)}\otimes\mathbb{S}_A^{(3,4)}\otimes\mathbb{S}_B^{(2,3)}\otimes\mathbb{S}_B^{(4,1)}$, then the fourth term can be written as
\begin{equation}
\begin{split}
\end{split}
\end{equation}
To write these six terms into one, we need to define some permutation operators, $\hat{W}_1=\mathbb{I}^{\otimes 4}$, $\hat{W}_2=\mathbb{I}^{\otimes 2}\otimes\mathbb{S}^{(3,4)}$, $\hat{W}_3=\mathbb{I}^{\otimes 2}\otimes\mathbb{S}^{(2,3)}$, $\hat{W}_4=\mathbb{I}^{\otimes 2}\otimes\mathbb{S}^{(2,4)}$, $\hat{W}_5=\mathbb{I}\otimes\overrightarrow{\Pi}^{(2,3,4)}$ and $\hat{W}_6=\mathbb{I}\otimes\overleftarrow{\Pi}^{(2,3,4)}$. Hence
\begin{equation}
\begin{split}
\end{split}
\end{equation}
where $\hat{O}_4$ can be graphically represented as
\begin{eqnarray}
\hat{O}_4=
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_1.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_2.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_3.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_4.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_5.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_6.png}
\end{tabular},
\end{eqnarray}
where the cross sign is the SWAP operator, and different colors stands for different parties, then
\begin{equation}
\begin{split}
\tr(\hat{O}_4^2)=\sum_{i,j=1}^6\tr[\hat{W}_i^\dagger\hat{O}\hat{W}_i\hat{W}_j^\dagger\hat{O}\hat{W}_j].
\end{split}
\end{equation}
We summarize the calculation in Table \ref{table:O_4calculation}. According to the mathematical property of operator SWAP and identity, $\tr(\mathbb{S})=d$ and $\tr(\mathbb{I}^{\otimes 2})=d^2$, it is easy to find that
\begin{eqnarray}\label{eq:O_4square}
\tr[\hat{O}_4^2]\leq 36\tr[\hat{O}^2]=
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_1.png}
\end{tabular}
\right)\left(
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_1.png}
\end{tabular}
\right)\right]
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_4_01.png}
\end{tabular}
\right)
\end{eqnarray}

\begin{table}
\centering
\setlength{\tabcolsep}{4mm}{
\begin{tabular}{c||c|c|c|c|c|c}
\hline\hline
\rule{0pt}{28pt}
& \includegraphics[scale=0.08]{O_4_1} & \includegraphics[scale=0.08]{O_4_2} &\includegraphics[scale=0.08]{O_4_3}&\includegraphics[scale=0.08]{O_4_4}&\includegraphics[scale=0.08]{O_4_5}&\includegraphics[scale=0.08]{O_4_6} \\
\hline\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.08]{O_4_1}&\includegraphics[scale=0.08]{O_4_01} &\includegraphics[scale=0.08]{O_4_02} &\includegraphics[scale=0.08]{O_4_03} &\includegraphics[scale=0.08]{O_4_04} &\includegraphics[scale=0.08]{O_4_05} &\includegraphics[scale=0.08]{O_4_06} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.08]{O_4_2}&\includegraphics[scale=0.08]{O_4_02} &\includegraphics[scale=0.08]{O_4_01} &\includegraphics[scale=0.08]{O_4_11} &\includegraphics[scale=0.08]{O_4_12} &\includegraphics[scale=0.08]{O_4_13} &\includegraphics[scale=0.08]{O_4_14} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.08]{O_4_3}&\includegraphics[scale=0.08]{O_4_03} &\includegraphics[scale=0.08]{O_4_11} &\includegraphics[scale=0.08]{O_4_01} &\includegraphics[scale=0.08]{O_4_21} &\includegraphics[scale=0.08]{O_4_22} &\includegraphics[scale=0.08]{O_4_23} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.08]{O_4_4}&\includegraphics[scale=0.08]{O_4_04} &\includegraphics[scale=0.08]{O_4_12} &\includegraphics[scale=0.08]{O_4_21} &\includegraphics[scale=0.08]{O_4_01} &\includegraphics[scale=0.08]{O_4_31} &\includegraphics[scale=0.08]{O_4_32} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.08]{O_4_5}&\includegraphics[scale=0.08]{O_4_05} &\includegraphics[scale=0.08]{O_4_13} &\includegraphics[scale=0.08]{O_4_22} &\includegraphics[scale=0.08]{O_4_31} &\includegraphics[scale=0.08]{O_4_01} &\includegraphics[scale=0.08]{O_4_41} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.08]{O_4_6}&\includegraphics[scale=0.08]{O_4_06} &\includegraphics[scale=0.08]{O_4_14} &\includegraphics[scale=0.08]{O_4_23} &\includegraphics[scale=0.08]{O_4_32} &\includegraphics[scale=0.08]{O_4_41} &\includegraphics[scale=0.08]{O_4_01} \\
\hline\hline
\end{tabular}}
\caption{ Multiplication of the terms in $\hat{O}_4^2$. In this table, we use blue and red colors to represent the two parties. ``x'' form signs stand for SWAP operators, and the horizontal lines stand for identity operators.  }
\label{table:O_4calculation}
\end{table}

The other three operators $\hat{O}_1$, $\hat{O}_2$ and $\hat{O}_3$ can be constructed using $\hat{O}_4$. Start from the third term,
\begin{equation}
\begin{split}
=&\tr[\hat{O}_4(\rho\otimes\mathbb{I}^{\otimes 3})(\mathbb{I}\otimes\hat{\rho}_1\otimes\hat{\rho}_2\otimes\hat{\rho}_3)]\\
=&\tr\left\{\tr_1[\hat{O}_4(\rho\otimes\mathbb{I}^{\otimes 3})](\hat{\rho}_1\otimes\hat{\rho}_2\otimes\hat{\rho}_3)\right\}\\
\end{split}
\end{equation}
so that
\begin{equation}
\begin{split}
\hat{O}_3=\tr_1[\hat{O}_4(\rho\otimes\mathbb{I}^{\otimes 3})].
\end{split}
\end{equation}
In order to make the graphic representation more convenient, when analyzing $\hat{O}_3$ $\hat{O}_2$ and $\hat{O}_1$, we will treat $\rho$ \textbf{as if} it is the tensor product of two reduced density matrices, $\rho=\rho_A\otimes\rho_B$, so that we can use two separated boxes to represent it. But one needs to remember that it is actually not the case. Therefore,
\begin{eqnarray}
\hat{O}_3&=\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_1_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_2_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_3_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_4_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_5_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_6_.png}
\end{tabular}
\\
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_1.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_2.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_3.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_4.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_5.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_6.png}
\end{tabular},
\end{eqnarray}
where we use the colored boxes to represent $\rho_A$ and $\rho_B$. To find out $\tr(\hat{O}_3^2)$, we make a similar table, Table.~\ref{table:O_3calculation}. By definition, taking trace of the terms in this table gives
\begin{eqnarray}
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_a1.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_a2.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_a3.png}
\end{tabular}
\right)=d_A^2\tr(\rho_A^2)
\end{eqnarray}
\begin{eqnarray}
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_b1.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_b2.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_b3.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_b4.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_b5.png}
\end{tabular}
\right)
\tr\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_3_b6.png}
\end{tabular}
\right)
\end{eqnarray}
\begin{eqnarray}\label{eq:O_3square}
\tr[\hat{O}_3^2]\leq 36\tr\left[
\left(
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_3_1.png}
\end{tabular}
\right)
\left(
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_3_1.png}
\end{tabular}
\right)
\right]
\begin{tabular}{c}
\includegraphics[scale=0.08]{O_3_11.png}
\end{tabular}
\right)
\end{eqnarray}

\begin{table}
\centering
\setlength{\tabcolsep}{3mm}{
\begin{tabular}{c||c|c|c|c|c|c}
\hline\hline
\rule{0pt}{28pt}
& \includegraphics[scale=0.1]{O_3_1} & \includegraphics[scale=0.1]{O_3_2} &\includegraphics[scale=0.1]{O_3_3}&\includegraphics[scale=0.1]{O_3_4}&\includegraphics[scale=0.1]{O_3_5}&\includegraphics[scale=0.1]{O_3_6} \\
\hline\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.1]{O_3_1}&\includegraphics[scale=0.1]{O_3_11} &\includegraphics[scale=0.1]{O_3_12} &\includegraphics[scale=0.1]{O_3_13} &\includegraphics[scale=0.1]{O_3_14} &\includegraphics[scale=0.1]{O_3_15} &\includegraphics[scale=0.1]{O_3_16} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.1]{O_3_2}&\includegraphics[scale=0.1]{O_3_21} &\includegraphics[scale=0.1]{O_3_22} &\includegraphics[scale=0.1]{O_3_23} &\includegraphics[scale=0.1]{O_3_24} &\includegraphics[scale=0.1]{O_3_25} &\includegraphics[scale=0.1]{O_3_26} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.1]{O_3_3}&\includegraphics[scale=0.1]{O_3_31} &\includegraphics[scale=0.1]{O_3_32} &\includegraphics[scale=0.1]{O_3_33} &\includegraphics[scale=0.1]{O_3_34} &\includegraphics[scale=0.1]{O_3_35} &\includegraphics[scale=0.1]{O_3_36} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.1]{O_3_4}&\includegraphics[scale=0.1]{O_3_41} &\includegraphics[scale=0.1]{O_3_42} &\includegraphics[scale=0.1]{O_3_43} &\includegraphics[scale=0.1]{O_3_44} &\includegraphics[scale=0.1]{O_3_45} &\includegraphics[scale=0.1]{O_3_46} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.1]{O_3_5}&\includegraphics[scale=0.1]{O_3_51} &\includegraphics[scale=0.1]{O_3_52} &\includegraphics[scale=0.1]{O_3_53} &\includegraphics[scale=0.1]{O_3_54} &\includegraphics[scale=0.1]{O_3_55} &\includegraphics[scale=0.1]{O_3_56} \\
\hline
\rule{0pt}{28pt}
\includegraphics[scale=0.1]{O_3_6}&\includegraphics[scale=0.1]{O_3_61} &\includegraphics[scale=0.1]{O_3_62} &\includegraphics[scale=0.1]{O_3_63} &\includegraphics[scale=0.1]{O_3_64} &\includegraphics[scale=0.1]{O_3_65} &\includegraphics[scale=0.1]{O_3_66} \\
\hline\hline
\end{tabular}}
\caption{ Multiplication of the terms in $\hat{O}_3^2$. Here, colored boxes stand for density matrices corresponding to different subsystems. ``X'' form operators and the horizontal lines follow the same definition as Table.~\ref{table:O_4calculation}, others are tri-partite cyclic permutation operators. }
\label{table:O_3calculation}
\end{table}

$\hat{O}_2$ can be similarly constructed from $\hat{O}_4$:
\begin{equation}
\begin{split}
\hat{O}_2=\frac{1}{2}\tr_{1,2}\left[\hat{O}_4(\rho\otimes\rho\otimes \mathbb{I}^{\otimes2})\right]=\frac{1}{2}\tr_1\left[\hat{O}_3(\rho\otimes\mathbb{I}^{\otimes2})\right],
\end{split}
\end{equation}
which can be graphically demonstrated as
\begin{eqnarray}
\hat{O}_2
\frac{1}{2}\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_1_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_2_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_3_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_4_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_5_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_6_.png}
\end{tabular}
\right)\\
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_1.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_2.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_3.png}.
\end{tabular}
\end{eqnarray}
After the calculation of $\tr[\hat{O}_4^2]$ and $\tr[\hat{O}_3^2]$, one could realize that, the largest term contains most horizontal lines (identity operator $\mathbb{I}$) without box, which contributes $d$ when taking trace. As a result,
\begin{eqnarray}\label{eq:O_2square}
\tr[\hat{O}_2^2]\leq 9\max\left\{\tr\left[
\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_1.png}
\end{tabular}
\right)^2
\right],
\tr\left[
\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_2_2.png}
\end{tabular}
\right)^2
\right]
\right\}=\max\left\{d_A^2M_4^{(2,3)},d_B^2M_4^{(2,3)}\right\}.
\end{eqnarray}

\begin{equation}
\begin{split}
\hat{O}_1=\frac{1}{6}\tr_{1,2,3}\left[\hat{O}_4(\rho^{\otimes 3}\otimes\mathbb{I})\right]=\frac{1}{6}\tr_{1,2}\left[\hat{O}_3(\rho^{\otimes 2}\otimes\mathbb{I})\right]=\frac{1}{3}\tr_1\left[\hat{O}_2(\rho\otimes\mathbb{I})\right],
\end{split}
\end{equation}
which can be graphically represented as
\begin{eqnarray}
\hat{O}_1 =
\frac{1}{3}\left(
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_1_1_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_1_2_.png}
\end{tabular}
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_1_3_.png}
\end{tabular}
\right)
\begin{tabular}{c}
\includegraphics[scale=0.1]{O_1.png}.
\end{tabular}
\end{eqnarray}
Notice that all lines in $\hat{O}_1$ are connected with boxes, so it is easily proved that
\begin{equation}\label{eq:O_1square}
\begin{split}
\tr[\hat{O}_1^2]\le 1.
\end{split}
\end{equation}

Now, substituting Eq.~\eqref{eq:O_4square}, Eq.~\eqref{eq:O_3square}, Eq.\eqref{eq:O_2square}, Eq.~\eqref{eq:O_1square}, and Lemma.~\ref{fact:shadow variance} into Eq.~\eqref{eq:M_4variance}, one can get the upper bound of the variance of the global shadow protocol
\begin{equation}\label{eq:GS variance}
\begin{split}
\mathrm{Var}\left(\hat{M}_4^{(2,3)}\right)\le C_1\frac{1}{M}+C_2\frac{\max\{d_A^2,d_B^2\}M_4^{(2,3)}}{M^2}+C_3\frac{D^2\tr(\rho^2)}{M^3}+C_4\frac{D^4}{M^4}.
\end{split}
\end{equation}
$C_1$, $C_2$, $C_3$ and $C_4$ in this equation are constants independent with number of snapshots $M$ and system dimension $D$.  They are actually not the same constants as those in Eq.~\eqref{eq:M_4variance}, but for simplicity, we use the same notation to represent them. To calculate the variance of local shadow protocol, we remark that
\begin{equation}
\begin{split}
\mathrm{locality}(\hat{O}_4)=4(N_A+N_B),\mathrm{locality}(\hat{O}_3)=3(N_A+N_B),\\
\mathrm{locality}(\hat{O}_2)=2(N_A+N_B),\mathrm{locality}(\hat{O}_1)=N_A+N_B.
\end{split}
\end{equation}
So we just need to make a little adjustment in the variance of global protocol to get the local one
\begin{equation}\label{eq:LS variance}
\begin{split}
\mathrm{Var}\left(\hat{M}_4^{(2,3)}\right)\le C_1\frac{D}{M}+C_2\frac{\max\{d_A^2,d_B^2\}M_4^{(2,3)}D^2}{M^2}+C_3\frac{D^5\tr(\rho^2)}{M^3}+C_4\frac{D^8}{M^4}.
\end{split}
\end{equation}
Therefore, in the large dimension scenario, the variance of global shadow protocol and the local one scale like $D^4/M^4$ and $D^8/M^4$ respectively. substituting this conclusion into the Chebyshev's equation, this proposition is proved. Similar analysis has also been conducted in the estimation of quantum negativity [\citenamefont{Elben et~al.}, 2020] and quantum fisher information [\citenamefont{Rath et~al.}, 2021].

\end{proof}

\subsection{Randomized Measurements Protocol and Statistical Analysis}\label{subsec:randomprotocol}
The logic of randomized measurement protocol is different with shadow estimation. We first need to write $M_2n^(2,3)(ρ_AB)$ in the observable form
\begin{equation}\label{eq:CCNR observable}
\begin{split}
M_{2n}^{(2,3)}(\rho_{AB})=\tr\left[(O_A\otimes O_B)\rho_{AB}^{\otimes 2n}\right],
\end{split}
\end{equation}
where $O_A=S_A^(1,2)⊗⋯⊗S_A^(2n-1,2n)$ and $O_B=S_B^(2n,1)⊗⋯⊗S_B^(2n-2,2n-1)$. Notice that the observable is just comprised by SWAP operators. SWAP operator, as the lowest order permutation operators, can be generated by averaging over 2-design global or qubit level local random unitary ensembles among virtual copies. Following the standard procedure of randomized measurements, we can design a protocol to measure $M_2n^(2,3)(ρ_AB)$:

\begin{algorithm}[H]
\caption{Measurement protocol of $M_{2n}^{(2,3)}$}
\label{algo:CCNR measurement protocol}
\begin{algorithmic}[1]
\Require
Sequentially prepared $2n\times N_U\times N_M$ $\rho_{AB}$
\Ensure
Probability distribution of the measurement outcomes conditioned on the evolution unitary $\mathrm{Pr}(\vec{s}_A,\vec{s}_B|U_A,U_B)$.

\For{$i= 1~\text{\textbf{to}}~N_U$} %\Comment{Enumerating all the locations}
\State Construct $2n$ unitary matrix $\{U_1=U_A^1\otimes U_B^1,U_2=U_A^1\otimes U_B^2,\cdots,U_{2n-1}=U_A^n\otimes U_B^n,U_{2n}=U_A^n\otimes U_B^1\}$ with each $U_A^i$ and $U_B^i$ sampled uniformly and independently from $N_A$ and $N_B$ qubits Clifford group (or constructed by tensor product of local unitaries $U_A^i=\bigotimes_{j=1}^{N_A}u_{j}$,$U_B^i=\bigotimes_{j=1}^{N_B}u_j$ with each $u_j$ sampled uniformly from qubit clifford group)
\State Operate these $2n$ unitaries on $\rho_{AB}$ to get the evolved states $\{U_1\rho_{AB}U_1^\dagger,\cdots,U_{2n}\rho_{AB}U_{2n}^\dagger\}$.
\For{$j= 1~\text{\textbf{to}}~N_M$}
\State  Measure the evolved states in the computational basis $\{\ket{\vec{s}_A,\vec{s}_B}\}$.
\State Record the measurement results.
\EndFor
\State Estimate the probabilities $\mathrm{Pr}(\vec{s}_A,\vec{s}_B|U_i)$.
\State Construct an unbiased estimator of $M_{2n}^{(2,3)}$ by these probabilities according to Proposition \ref{prop:CCNR measurement}

\EndFor
\State Take an average of the estimators in each round to get the final estimator of $M_{2n}^{(2,3)}$.
\end{algorithmic}
\end{algorithm}

\begin{proposition}\label{prop:CCNR measurement}
The unbiased estimator of $2n$-th moment of realignment matrix can be constructed using these probabilities $\mathrm{Pr}(\vec{s}_A,\vec{s}_B|U_i)$
\begin{equation}\label{eq:CCNR M2n}
\begin{split}
\hat{M}_{2n}^{(2,3)}&=\sum_{\substack{\vec{s}_A^1,\cdots,\vec{s}_A^{2n}\\ \vec{s}_B^1,\cdots,\vec{s}_B^{2n}}}\Omega
\prod_{i=1}^n X_A(\vec{s}_A^{2i-1},\vec{s}_A^{2i})X_B(\vec{s}_B^{2i},\vec{s}_B^{2i+1}),
\end{split}
\end{equation}
\begin{equation}\label{eq:CCNR omega}
\begin{split}
\Omega=\prod_{i=1}^n\mathrm{Pr}(\vec{s}_A^{2i-1},\vec{s}_B^{2i-1}|U_{2i-1})\mathrm{Pr}(\vec{s}_A^{2i},\vec{s}_B^{2i}|U_{2i})
\end{split}
\end{equation}
is the  product of probabilities, $U_i$ is chosen from the prepared unitary evolution group $\{U_1=U_A^1\otimes U_B^1,U_2=U_A^1\otimes U_B^2,\cdots,U_{2n-1}=U_A^n\otimes U_B^n,U_{2n}=U_A^n\otimes U_B^1\}$ as stated in Algo.~\ref{algo:CCNR measurement protocol}, $X(\vec{s},\vec{s}')$ is the weight function. For globally constructed $U_{Ai}$ and $U_{Bi}$
\begin{equation}
\begin{split}
\end{split}
\end{equation}
for locally constructed $U_A^i=\bigotimes_{j=1}^{N_A}u_{j}$,$U_B^i=\bigotimes_{j=1}^{N_B}u_j$
\begin{equation}
\begin{split}
\end{split}
\end{equation}
where $|\vec{s}|$ is the length of $\vec{s}$, $D[\vec{s},\vec{s}']$ is the Hamming distance between $\vec{s}$ and $\vec{s}'$. In Eq.~\eqref{eq:CCNR M2n}, we set $\vec{s}_B^{2n+1}=\vec{s}_B^1$.

The estimator constructed above satisfies
\begin{equation}
\begin{split}
\mathbb{E}_{U,\vec{s}}\left(\hat{M}_{2n}^{(2,3)}\right)&=M_{2n}^{(2,3)}
\end{split}
\end{equation}
where $\mbb{E}_{U,\vec{s}}$ denotes the expectation over all Clifford groups and measurement results.

\end{proposition}

\begin{proof}
To prove this proposition, we define the postprocessing operator $X_{g/l}=\sum_{\vec{s},\vec{s}'}X_{g/l}(\vec{s},\vec{s}')\ketbra{\vec{s},\vec{s}'}{\vec{s},\vec{s}'}$, which can be used to generate SWAP operator [\citenamefont{Brydges et~al.}, 2019, \citenamefont{Elben et~al.}, 2019] through a twirling channel,
\begin{equation}\label{eq:SWAP construct global}
\begin{split}
\Phi_2(X_g)=\mathbb{E}_{U\in\mathcal{E}}U^{\otimes 2}X_g U^{\dagger\otimes 2}=\mathbb{S},
\end{split}
\end{equation}
where $\mathcal{E}$ is a unitary 2-design, which is fulfilled by Clifford group.  Substituting Born's rule, $\mathbb{E}_{\vec{s}}\mathrm{Pr}(\vec{s}_A,\vec{s}_B|U_A,U_B)=\tr\left[\ketbra{\vec{s}_A,\vec{s}_B}{\vec{s}_A,\vec{s}_B}(U_A\otimes U_B)\rho_{AB}(U_A^\dagger\otimes U_B^\dagger)\right]$, into the right hand side of Eq.~\eqref{eq:CCNR M2n}, and consider the global case,
\begin{equation}
\begin{split}
R.H.S&=\sum_{\substack{\vec{s}_A^1,\cdots,\vec{s}_A^{2n}\\ \vec{s}_B^1,\cdots,\vec{s}_B^{2n}}}\mbb{E}_{U}\prod_{i=1}^n X_A(\vec{s}_A^{2i-1},\vec{s}_A^{2i})X_B(\vec{s}_B^{2i},\vec{s}_B^{2i+1})\\
&\tr\left\{\left(\bigotimes_{i=1}^n\ketbra{\vec{s}_A^{2i-1}\vec{s}_A^{2i}}{\vec{s}_A^{2i-1}\vec{s}_A^{2i}}\otimes\ketbra{\vec{s}_B^{2i-1}\vec{s}_B^{2i}}{\vec{s}_B^{2i-1}\vec{s}_B^{2i}}\right)\left(\bigotimes_{i=1}^n U_i\right)\rho_{AB}^{\otimes 2n}\left(\bigotimes_{i=1}^n U_i\right)^\dagger\right\}\\
&=\mbb{E}_{U}\tr\left\{\left(\bigotimes_{i=1}^n X_A^{(2i-1,2i)}\otimes X_B^{(2i,2i+1)}\right)\left(\bigotimes_{i=1}^n U_i\right)\rho_{AB}^{\otimes 2n}\left(\bigotimes_{i=1}^n U_i\right)^\dagger\right\}\\
&=\mbb{E}_{U}\tr\left\{\left(\bigotimes_{i=1}^n U_i\right)^\dagger\left(\bigotimes_{i=1}^n X_A^{(2i-1,2i)}\otimes X_B^{(2i,2i+1)}\right)\left(\bigotimes_{i=1}^n U_i\right)\rho_{AB}^{\otimes 2n}\right\}\\
&=\tr\left\{\left(\bigotimes_{i=1}^n \Phi_2(X_A^{(2i-1,2i)})\otimes \Phi_2(X_B^{(2i,2i+1)})\right)\rho_{AB}^{\otimes 2n}\right\}\\
&=\tr\left\{\left(\bigotimes_{i=1}^n \mathbb{S}_A^{(2i-1,2i)}\otimes\mathbb{S}_B^{(2i,2i+1)} \right)\rho_{AB}^{\otimes 2n}\right\}\\
\end{split}
\end{equation}
Here we use subscript to denote the party and superscript to denote the copy number, like $X_A^{(2i-1,2i)}$ is the operator acting on the $A$ parties of the $(2i-1)$-th and $2i$-th copies of $\rho_{AB}$.

The proof of local case is quite similar to global one. Notice that $X_l=\bigotimes_{i=1}^N X_{gi}$, hence
\begin{equation}\label{eq:SWAP construct local}
\begin{split}
\mathbb{E}_U\left(\bigotimes_{i=1}^Nu_i\right)^{\otimes 2}X_l\left(\bigotimes_{i=1}^Nu_i\right)^{\dagger\otimes 2}=\mathbb{E}_U\left(\bigotimes_{i=1}^Nu_i^{\otimes 2}X_{gi}u_i^{\dagger \otimes 2}\right)=\bigotimes_{i=1}^N\Phi_2(X_{gi})=\bigotimes_{i=1}^N\mathbb{S}_i=\mathbb{S}.
\end{split}
\end{equation}
Then, Eq.~\eqref{eq:CCNR M2n} also holds for local unitary.

\end{proof}

\begin{figure}[htbp]
\centering