# Auxiliary-free replica shadow estimation

Qing Liu  0000-0002-1576-1975 These authors contributed equally to this work.
Key Laboratory for Information Science of Electromagnetic Waves (Ministry of
Education), Fudan University, Shanghai 200433, China  Zihao Li
0000-0002-2766-7521 These authors contributed equally to this work. State Key
Laboratory of Surface Physics and Department of Physics, Fudan University,
Shanghai, 200433, China  Institute for Nanoelectronic Devices and Quantum
Computing, Fudan University, Shanghai, 200433, China  Center for Field Theory
and Particle Physics, Fudan University, Shanghai, 200433, China  Xiao Yuan
0000-0003-0205-6545 Center on Frontiers of Computing Studies, Peking
University, Beijing 100871, China  School of Computer Science, Peking
University, Beijing 100871, China  Huangjun Zhu  0000-0001-7257-0764 State Key
Laboratory of Surface Physics and Department of Physics, Fudan University,
Shanghai, 200433, China  Institute for Nanoelectronic Devices and Quantum
Computing, Fudan University, Shanghai, 200433, China  Center for Field Theory
and Particle Physics, Fudan University, Shanghai, 200433, China  You Zhou
0000-0003-0886-077X<EMAIL_ADDRESS>Key Laboratory for Information
Science of Electromagnetic Waves (Ministry of Education), Fudan University,
Shanghai 200433, China

###### Abstract

Efficiently measuring nonlinear properties, like the entanglement spectrum, is
a significant yet challenging task from quantum information processing to
many-body physics. Current methodologies often suffer from an exponential
scaling of the sampling cost or require auxiliary qubits and deep quantum
circuits. To address these limitations, we propose an efficient auxiliary-free
replica shadow (AFRS) framework, which leverages the power of the joint
entangling operation on a few input replicas while integrating the mindset of
shadow estimation. We rigorously prove that AFRS can offer exponential
improvements in estimation accuracy compared with the conventional shadow
method, and facilitate the simultaneous estimation of various nonlinear
properties, unlike the destructive swap test. Additionally, we introduce an
advanced local-AFRS variant tailored to estimating local observables with even
constant-depth local quantum circuits, which significantly simplifies the
experimental realization compared with the general swap test. Our work paves
the way for the application of AFRS on near-term quantum hardware, opening new
avenues for efficient and practical quantum measurements.

## Introduction

Quantum measurement is fundamental to quantum mechanics and plays a
significant role in the development of quantum technologies. Due to the
intrinsic linearity of quantum mechanics, the measurement of linear
observables, such as correlation functions and magnetization, is more
attainable compared with challenging nonlinear ones like purity, entanglement
entropy [1, 2, 3], and higher-order moments [4, 5, 6, 7]. These nonlinear
measurements have valuable applications in quantum foundations [8, 9], quantum
many-body physics [10, 11, 12, 13], quantum error mitigation, and algorithm
design [14, 15, 16].

There are mainly two methods to realize nonlinear measurements. The first is
the recently developed randomized measurement approach [17], particularly
shadow estimation [18, 19], which can return an unbiased estimator
$\hat{\rho}$ by repeating the projective measurement on the single-copy state
$\rho$ under randomly selected bases. Although the estimation of the nonlinear
function can be classically post-processed from many such single-shot
snapshots, this approach introduces heavy sampling and post-processing budgets
that scale exponentially with the system size [6, 7]. The second method,
called the (generalized) swap test, poses higher requirements on hardware. It
executes entangling measurement among a few copies of $\rho$ simultaneously
[1, 4], similar to the famous Bell measurement on two replicas. This method is
much more efficient regarding sampling time and post-processing than the first
one for measuring nonlinear properties [20, 21].

Figure 1:  An illustration of the auxiliary-free replica shadow estimation
framework. In the experimental stage, one prepares $t$ replicas $\rho^{\otimes
t}$, and evolves each $\rho$ with the same $V\in\mathcal{E}$, where
$\mathcal{E}$ can be global random Clifford ensemble
$\mathcal{E}_{\mathrm{Cl}}$ or tensor-product of single-qubit ones
$\mathcal{E}_{\mathrm{LCl}}$. After the joint unitary $\mathcal{R}$, one
measures all replicas in the computational basis, which returns the outcome
$\mathbf{x}$. During the classical post-processing stage, the outcome
$\mathbf{x}$ undergoes a random mapping to $\mathbf{b}$ via the strategy
$\mathcal{G}_{\mathbf{x}\to\mathbf{b}}$, while concurrently accessing the
function $f(\mathbf{x})$ as defined in 1. By Eq. 4 in Theorem 1, one can
construct the unbiased estimator of $\rho^{t}$ as $\widehat{\rho^{t}}$.
Repeating the whole procedure for $M$ times, one can obtain the shadow set
$\\{\widehat{\rho^{t}_{(i)}}\\}_{i=1}^{M}$, which would lead to diverse
applications from the quantum foundation and information processing to many-
body physics, such as estimating the nonlinear properties $\tr(O\rho^{t})$ of
various observables in a multiplex way for virtual distillation.

The swap test in general applies a (generalized) controlled-swap circuit
before the final native measurements, necessitating one auxiliary qubit and a
circuit depth that scales with the qubit number $n$. In certain cases, we can
reduce the circuit depth to a constant, but this would inevitably require
$\mathcal{O}(n)$ auxiliaries initialized in an entangled state [22]. These
requirements make its realization on near-term quantum platforms challenging,
and existing works are limited to a few small-scale demonstrations [23, 24,
25, 26, 27]. Alternatively, entangling measurements without any auxiliaries,
such as the destructive swap test or its variants [14, 15, 28, 29], have been
realized in cold-atom [2, 10, 14, 30] and superconducting-qubit systems [20,
31]. Nevertheless, it is a widely held belief that the information acquired
from such destructive measurements is highly limited compared with the full
swap test [15, 32, 21]. For instance, consider the $t$-th order nonlinear
function $\tr(O\rho^{t})$, which is essential in virtual distillation [14] and
quantum error mitigation tasks [15, 16]. Current state-of-the-art developments
are restricted to the moment function $\tr(\rho^{t})$ (say, $O=\mathbb{I}$) or
cases where the observable is some single-qubit ones [15] with case-by-case
tailored circuits. Realizing general measurements of complex nonlinear
functions like $\tr(O\rho^{t})$ without the expensive swap test remains an
interesting open problem.

Here, we address the problem by proposing a general framework called
auxiliary-free replica shadow (AFRS) estimation, as illustrated in Fig. 1.
AFRS inherits the mindset of shadow estimation and elaborately integrates it
with fine-tuned entangling measurements and efficient classical post-
processings. In particular, it first applies the _same_ random unitary
evolution for all $t$ replicas, and then conducts an entangling measurement
among them, which is realized by a joint unitary $\mathcal{R}$ followed by
native measurements. Through a delicately designed classical (random) mapping
strategy, AFRS returns the estimator $\widehat{\rho^{t}}$ along the shadow
estimation procedure, which can be used as the estimator
$\tr\\!\big{(}O\widehat{\rho^{t}}\big{)}$ for different $O$, and other more
complicated nonlinear functions. Notably, the variance or statistical
estimation error of AFRS shows exponential advantages compared with the
original (single-copy) shadow method [19, 33], and could even become constant
and independent of the system size for certain cases. By simplifying the
circuit compilation of $\mathcal{R}$, we reduce the circuit depth to a
constant in important scenarios of estimating general local observables, such
as many-body Hamiltonians, enabling practical applications on near-term
platforms. Moreover, thanks to the shadow nature of the framework, the
randomized measurement data can be multiplexed, making the sampling time
scales logarithmically $\mathcal{O}(\log(L))$ with the total number $L$ of
observables, i.e., an exponential advantage compared with previous methods
[15, 16]. We anticipate the proposed AFRS, an efficient method for estimating
nonlinear functions of quantum objects, would serve as a basic and useful
subroutine in many quantum information and computing tasks.

## Results

We begin by illustrating the general protocol of AFRS and its performance
guarantee in Section .1. Different from the usually-used Bell measurements,
the general form of the essential entangling unitary $\mathcal{R}$ among
replicas are given in Section .2. We further show the explicit circuit
compilation of $\mathcal{R}$ in the general qudit and many-qubit systems in
Section .3, with ideas and techniques stimulated from quantum error correction
code [34], which may be of independent interest. In particular, $\mathcal{R}$
can be realized with controlled-NOT and Hadamard gates. We develop the local-
AFRS protocol in Section .4, tailoring to estimating local observables with a
significantly reduced circuit depth. In Section .5, we demonstrate practical
advantages of local-AFRS in the virtual distillation task.

###  .1 Framework of AFRS

Let us first give a brief introduction to the original shadow (OS) estimation
[19]. We consider a general quantum state $\rho\in\mathcal{H}_{d}$, and use
$\ket{\mathbf{b}}$ with $\mathbf{b}\in[d]$ to denote the corresponding
computational basis. If it is an $n$-qubit system, say
$\mathcal{H}_{d}=\mathcal{H}_{2}^{\otimes n}$, we expand the basis to the
binary form $\mathbf{b}=\\{b_{1},b_{2}\cdots,b_{n}\\}$ for every qubit. OS
evolves the quantum state $\rho$ with some unitary sampled from some random
ensemble $V\in\mathcal{E}$, and measures it in the computational basis to get
$\mathbf{b}$ with the Born probability
$\Pr(\mathbf{b}|V)=\bra{\mathbf{b}}V(\rho)\ket{\mathbf{b}}$. Here,
$V(\rho):=V\rho V^{{\dagger}}$ for short. One can then construct a one-shot
unbiased estimator
$\hat{\rho}=\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)$
on a classical computer. The classical inverse channel $\mathcal{M}^{-1}$
depends on the chosen unitary ensemble, with the $n$-qubit Clifford unitary
$\mathcal{E}_{\mathrm{Cl}}$ and tensor-product of single-qubit Clifford
$\mathcal{E}_{\mathrm{LCl}}$, as two typical choices [19, 35, 36, 37]. One
repeats such “measurement $+$ prepare” procedure for $M$ times to collect the
shadow set $\\{\hat{\rho}_{i}\\}_{i=1}^{M}$. The expectation value of $O$ can
be estimated via $\tr(O\hat{\rho})$, which can be efficient for local
observable and fidelity estimation. On the other hand, nonlinear functions
like $\tr(O\rho^{t})$ can be estimated by
$\tr(O\hat{\rho}_{1}\hat{\rho}_{2}\cdots\hat{\rho}_{t})$, with $t$ distinct
snapshots from the shadow set. One can reduce the estimation error by
increasing the size of the shadow set. However, for such kind of nonlinear
functions, even for $O$ being identity or some local Pauli operators, the
variance of the estimation or equivalently the sampling time would inevitably
become exponentially large with the qubit number [19, 6], which indicates a
prohibitive sampling time even for a medium-scale system. Entangling
operations on multiple replicas simultaneously are believed to enhance the
sampling performance in the estimation, even showing an exponential advantage,
such as the estimation of purity with the Bell measurement [38, 20].

Here, we propose the framework of AFRS, which borrows the mindset of
randomized measurements and shadow estimation, and tactfully incorporates them
with the entangling operation on all $t$ replicas, say $\rho^{\otimes t}$, to
faithfully return the unbiased estimator $\widehat{\rho^{t}}$ of the monomial
of the quantum state, as shown in Fig. 1. Note that other nonlinear functions
can be directly estimated using our framework, in particular, the $t$ input
replicas can be different, and can also be with some additional pre- or post-
operations [21]. Such nonlinear functions include the state overlap and
fidelity [39, 40], the purity and higher-order moments [41, 42], quantum
Fisher information [43, 44], out-of-time-ordered correlators [11, 13],
topological invariants [12, 45], etc. As a result, hereafter we mainly adopt
the monomial function $\tr(O\rho^{t})$ for the following discussion.

For the $t$-replica Hilbert space, $\mathcal{H}_{d}^{\otimes t}$, we use
$\\{|\mathbf{x}\rangle:\mathbf{x}=\\{\mathbf{x}_{1}\mathbf{x}_{2}\cdots\mathbf{x}_{t}\\}\in\mathbb{B}(d,t)\\}$
as the computational basis, with $\mathbb{B}(d,t):=[d]^{t}$. For the case of a
subsystem $A\in[n]$ of the $n$-qubit system, we use $\mathbf{b}^{A}$ labels
the bit-string restricted on $A$, and similarly
$\mathbf{x}^{A}=\\{\mathbf{x}_{1}^{A}\mathbf{x}_{2}^{A}\cdots\mathbf{x}_{t}^{A}\\}$
for the bit-strings restricted on $A$ for all replicas. For our problem, the
key point is to generate the corresponding conditional ‘probability’

$\Pr(\mathbf{b}|V)=\bra{\mathbf{b}}V(\rho^{t})\ket{\mathbf{b}},$ (1)

by substituting $\rho^{t}$ for $\rho$ as in OS.

Note that $\rho^{t}$ in general is not a legal quantum state, and thus Eq. 1
is not a real probability. However, we can use entangling operation to
generate such _fake_ probability by the following observation. Define the
observable in the computational basis as

$Q_{\mathbf{b}}=t^{-1}\sum_{i=1}^{t}\outerproduct{\mathbf{b}}{\mathbf{b}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)},$
(2)

which is symmetric among all $t$ replicas, and only takes $\mathbf{b}$ on one
of them. It is not hard to check that such a symmetric observable commutes
with the $t$-copy shift operation, i.e., $[S_{t},Q_{\mathbf{b}}]=0$. Here,
$S_{t}$ is the clockwise shift operation, with
$S_{t}\ket{\mathbf{x}_{1}\mathbf{x}_{2}\cdots\mathbf{x}_{t}}=\ket{\mathbf{x}_{2}\mathbf{x}_{3}\cdots\mathbf{x}_{1}}$.
As a result, $S_{t}$ and $Q_{\mathbf{b}}$ can be diagonalized simultaneously
in some common orthonormal basis, denoted by $\\{\ket{\psi_{\mathbf{x}}}\\}$.

######  Observation 1.

The fake conditional probability in Eq. 1 can be written as

$\displaystyle\Pr(\mathbf{b}|V)=\tr\left(S_{t}Q_{\mathbf{b}}V(\rho)^{\otimes
t}\right)$
$\displaystyle=\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}$
$\displaystyle=\sum_{\mathbf{x}}f(\mathbf{x})\Pr(\mathbf{b}|\mathbf{x})\Pr(\mathbf{x}|V).$
(3)

Here, $f(\mathbf{x}):=\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}}$,
which is a function depending the measurement result $\mathbf{x}$,
$\Pr(\mathbf{b}|\mathbf{x}):=\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}$
and $\Pr(\mathbf{x}|V):=\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}$, respectively, which are _legal_ conditional
probabilities.

The proof of 1 is left in Section A.1. It is not hard to check that
$\sum_{\mathbf{b}}\Pr(\mathbf{b}|\mathbf{x})=1$, and we call such conditional
probability as the mapping strategy,
$\mathcal{G}_{\mathbf{x}\to\mathbf{b}}:=\Pr(\mathbf{b}|\mathbf{x})$, since it
is generated by the classical random assignment, which is explicitly shown in
Section .2 where the form of $\\{\ket{\psi_{\mathbf{x}}}\\}$ is also given. On
the other hand, $\Pr(\mathbf{x}|V)$ can be generated by the quantum way, i.e.,
measuring $V(\rho)^{\otimes t}$ in the basis $\\{\ket{\psi_{\mathbf{x}}}\\}$.
In principle, one can find a unitary operation $\mathcal{R}$ on
$\mathcal{H}_{d}^{\otimes t}$ such that
$\mathcal{R}\ket{\psi_{\mathbf{x}}}=\ket{\mathbf{x}}$. In this way, the whole
quantum circuit is summarized in Fig. 1.

Now we are in the position to give our main protocol in Theorem 1, which aims
to conduct the classical shadow of $\rho^{t}$ from the randomized and
entangling measurement results from Fig. 1, and a rigorous performance
guarantee for the corresponding estimation variance.

######  Theorem 1.

Suppose one conducts the circuit proposed in Fig. 1 for once (i.e., single-
shot, $M=1$), the unbiased estimator of $\rho^{t}$ shows

$\widehat{\rho^{t}}:=\mathrm{Re}(f(\mathbf{x}))\cdot\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right),$
(4)

such that
$\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\widehat{\rho^{t}}=\rho^{t}$. Here,
the inverse channel $\mathcal{M}^{-1}$ depends on the random unitary ensemble
$\mathcal{E}$; $\mathbf{x}$ is the actual measurement outcome from the quantum
circuit in Fig. 1; the function $f$, and the mapping result $\mathbf{b}$ are
given by 1.

The unbiased estimator of some observable $O$ on $\rho^{t}$, $\tr(O\rho^{t})$,
is $\widehat{o_{t}}=\tr(O\widehat{\rho^{t}})$, and its variance shows

$\displaystyle\mathrm{Var}(\widehat{o_{t}})\leq\left\|O_{0}\right\|_{\mathrm{sh},\mathcal{E}}^{2}+\|O\|_{\infty}^{2},$
(5)

where $O_{0}=O-\operatorname{tr}(O)\mathbb{I}_{d}/d$ is the traceless part of
$O$; and $\|O_{0}\|_{\mathrm{sh},\mathcal{E}}$ is the shadow norm depending on
the unitary ensemble $\mathcal{E}$ [19].

The proof of Theorem 1 is left in Section A.2. Indeed, one can also construct
the unbiased estimator $\hat{\rho}$ of $\rho$ by just ignoring the function
term $\mathrm{Re}(f(\mathbf{x}))$ in Eq. 4. We remark that the shadow set
$\\{\widehat{\rho^{t}_{(i)}}\\}_{i=1}^{M}$ is collected by repeating the
single-shot procedure for $M$ times. Then the final variance to estimate
$o_{t}$ in Eq. 5 would be reduced by a factor of $M$ via averaging these
i.i.d. estimators. The shadow norm introduced in Ref. [19], is used here to
quantify the single-shot statistical variance as shown in Eq. 5, and thus a
larger variance leads to a large sampling time $M$ in the practical
estimation. Even though we show the result for one observable $O$ here, the
shadow set can be reused to simultaneously estimate a set of observables
$\\{O^{i}\\}$. By using the median-of-means technique as that in OS [19], the
sampling time $M$ or the size of the shadow set only scales logarithmically
with the total number of the observables, and we will discuss this advantage
in detail in Section .4.

Theorem 1 indicates the significant reduction in the variance of estimating
nonlinear functions using AFRS compared with the one using OS. In particular,
the (square of) shadow norm $\|O_{0}\|^{2}_{\mathrm{sh}}$ is used as an upper
bound of the variance for estimating linear observable $O$ with $\hat{\rho}$
in OS. Now by Eq. 5 and ignoring the second term $\|O\|_{\infty}^{2}$ which is
usually a constant, it is clear to see that AFRS reduces the variance for
estimating nonlinear function $\tr(O\rho^{t})$ to the linear case, which leads
to exponential advantages in many cases of interest. Take the random local
Clifford ensemble as example, where $V\in\mathcal{E}_{\mathrm{LCl}}$, the
shadow norm of a $\omega$-local observable $O$ is upper bounded by
$\left\|O_{0}\right\|_{\mathrm{sh},\mathrm{LCl}}^{2}\leq
4^{\omega}\|O\|_{\infty}^{2}$. Thus, according to Eq. 5, the variance of
estimating $\tr(O\rho^{t})$ is independent of the qubit number $n$ using AFRS;
However, for OS, as $\tr(O\rho^{t})=\tr(O_{t}\rho^{\otimes t})$ with
$O_{t}:=(O\otimes{\mathbb{I}}_{d}^{t-1})S_{t}$ being a highly global
observable due to $S_{t}$, the shadow norm of $O_{t}$ and thus the estimation
variance would scale exponentially with $n$, which is demonstrated numerically
in Fig. 2 (a). Similarly, for random Clifford ensemble
$V\in\mathcal{E}_{\mathrm{Cl}}$, the shadow norm is
$\left\|O_{0}\right\|_{\mathrm{sh},\mathrm{Cl}}^{2}\leq 3\tr(O^{2})$, while
$\left\|O_{t}\right\|_{\mathrm{sh},\mathrm{Cl}}^{2}\sim d^{t-1}\tr(O^{2})$,
which implies an exponential scaling advantage using the AFRS protocol for the
fidelity estimation with $O$ a pure state, as shown in Fig. 2 (b).

Figure 2: Scaling of estimation error with respect to the qubit number for OS,
AFRS and local-AFRS protocols. The quantum state is a noisy $n$-qubit GHZ
state with $\rho=0.7\ket{\text{GHZ}}\bra{\text{GHZ}}+0.3{\mathbb{I}}_{d}/d$.
In (a), $V$ is sampled from the local Clifford ensemble
$V\in\mathcal{E}_{\mathrm{LCl}}$ for estimating $o_{2}=\tr(O\rho^{2})$, with
$O=Z_{1}Z_{2}$ being a local Pauli observable. In (b), random Clifford
ensemble, where $V\in\mathcal{E}_{\mathrm{Cl}}$, is utilized for estimating
the fidelity $F_{2}=\bra{\text{GHZ}}\rho^{2}\ket{\text{GHZ}}$ with
$O=\ket{\text{GHZ}}\bra{\text{GHZ}}$. In (a), we take the repetition number
$M=50$ for the AFRS and local-AFRS protocols, and take $M=100$ for OS. While
in (b), we take $M=10$ for AFRS and $M=20$ for OS. Thus, the estimation errors
with different protocols are compared under the same sample number of $\rho$.

###  .2 The matrix form of $\mathcal{R}$

In this section, we give the matrix form of the unitary $\mathcal{R}$
expressed in the computational basis $\\{|\mathbf{x}\rangle\\}$. Two strings
$\mathbf{x},\mathbf{y}\in\mathbb{B}(d,t)$ are said to be in the same class if
there is an integer $r$ such that $\tau^{r}(\mathbf{x})=\mathbf{y}$. That is,
they can be transformed to each other by the one-step cyclic permutation
$\tau$ which corresponds to the unitary $S_{t}$ on $t$ copies, i.e.,
$S_{t}\ket{\mathbf{x}}=\ket{\tau(\mathbf{x})}$. For instance, the strings
$01353$ and $35301$ are in the same class of $\mathbb{B}(6,5)$. It is clear
that two distinct classes have no intersection, and $\mathbb{B}(d,t)$ is the
union of all distinct classes. For each class of $\mathbb{B}(d,t)$, we choose
its smallest element $\mathbf{z}$ as the representative, and use
$[\mathbf{z}]$ to represent the class. For example, $\mathbb{B}(2,4)$ can be
divided into six distinct classes, and one of them is
$[0011]=\\{0011,0110,1100,1001\\}$.

######  Lemma 1.

For any class $[\mathbf{z}]$ of $\mathbb{B}(d,t)$, we define its corresponding
normalized $t$-qudit states by Fourier transformation as

$\displaystyle\ket{\Psi_{[\mathbf{z}]}^{(k)}}=\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{rk\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\ket{\tau^{r}(\mathbf{z})}$ (6)

for $k=0,1,\dots,|[\mathbf{z}]|-1$, where $|[\mathbf{z}]|$ is the cardinality
of $[\mathbf{z}]$. The set of states
$\left\\{\ket{\Psi_{[\mathbf{z}]}^{(k)}}\right\\}_{[\mathbf{z}],k}$ forms an
orthonormal basis on $\mathcal{H}_{d}^{\otimes t}$. In addition,
$\ket{\Psi_{[\mathbf{z}]}^{(k)}}$ is the common eigenstate of $S_{t}$ and
$Q_{\mathbf{b}}$, with eigenvalue $\exp\big(\\!-\frac{k\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\big{missing})$ and $t^{-1}\times($the number of
$\mathbf{b}$ in the string $\mathbf{z})$, respectively.

Lemma 1 is proved in Appendix B. Based on the basis in Lemma 1, it is direct
to write down the corresponding entangling unitary as

$\mathcal{R}=\sum_{[\mathbf{z}]}\sum_{k=0,\dots,|[\mathbf{z}]|-1}\ket{\tau^{k}(\mathbf{z})}\bra{\Psi_{[\mathbf{z}]}^{(k)}}.$
(7)

In the case of $t=2$, the permutation $\tau$ of $S_{2}$ is just the swap
operation, and any string in $\mathbb{B}(d,t)$ can be written as $ab$ for
$a,b\in\\{0,\dots,d-1\\}$. If $a=b$, the class $[ab]=\\{ab\\}$; if $a<b$, the
class $[ab]=\\{ab,ba\\}$. In this way, the unitary $\mathcal{R}$ satisfies

$\displaystyle\mathcal{R}|ab\rangle=\begin{cases}|ab\rangle&a=b,\\\
\frac{1}{\sqrt{2}}(|ba\rangle+|ab\rangle)&a<b,\\\
\frac{1}{\sqrt{2}}(|ba\rangle-|ab\rangle)&a>b.\\\ \end{cases}$ (8)

and it is clear that $\mathcal{R}=\mathcal{R}^{{\dagger}}$ is Hermitian.

According to the eigenvalue of $S_{t}$ shown in Lemma 1, it is clear to see
that (the real part) of the function in 1 and Theorem 1 reads

$\mathrm{Re}(f(\mathbf{x}))=\cos(\frac{k\cdot 2\pi}{|[\mathbf{z}]|}),$ (9)

where $[\mathbf{z}]$ is the class that $\mathbf{x}$ belongs to, and the
integer $k$ satisfies $\tau^{k}(\mathbf{z})=\mathbf{x}$. For the case $t=2$,
we have $\mathrm{Re}(f(\mathbf{x}))=\pm 1$, which are the eigenvalues of
$S_{2}$.

It is also worth noting, from the eigenvalue of $Q_{\mathbf{b}}$ in Lemma 1,
that the mapping strategy $\mathcal{G}_{\mathbf{x}\to\mathbf{b}}$ is just to
randomly choose $\mathbf{b}$ from
$\\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{t}\\}$ of $\mathbf{x}$
with uniform probability $1/t$. In practise, instead of such a random
assignment, one can alternatively choose $\mathbf{b}$ as the average of all
$\mathbf{x}_{i}$. Such kind of treatment directly simulates the random
assignment process and inherits the multi-shot nature [46, 47], and would
further reduce the variance in Eq. 5 to a factor of $t^{-1}$ in some cases. We
give the form of the new estimator and the detailed variance analysis in
Appendix E.

###  .3 Circuit compilation of $\mathcal{R}$

After introducing the explicit matrix form of $\mathcal{R}$, here in this
section, we study the circuit compilation, which is significant for the
realization in quantum computing platforms. Here, we focus on the $t=2$
scenario, which is the most promising to demonstrate on near-term devices.
That is, $\mathcal{R}$ acting on $\mathcal{H}_{d}^{\otimes 2}$ is Hermitian
with matrix form in Eq. 8. We start from the single-qubit case with $d=2$, and
then extend to the general qudit case, and finally give the construction of
the most challenging $n$-qubit case with
$\mathcal{H}_{d}=\mathcal{H}_{2}^{\otimes n}$, which may be of independent
interest.

For the single-qubit case, as indicated by Eq. 8, $\mathcal{R}$ creates
superposition, say some Hadamard-like rotation, in the subspace
$\mathcal{V}={\rm span}(\ket{01},\ket{10})$, but is inoperative for the
complement subspace $\mathcal{V}^{\perp}={\rm span}(\ket{00},\ket{11})$. As
shown in Fig. 3 (a),

$\displaystyle\mathcal{R}=CX_{1\rightarrow 2}CH_{2\rightarrow
1}CX_{1\rightarrow 2}$ (10)

with $CX$ and $CH$ for the controlled-$X$ (say, CNOT) and controlled-$H$ gate,
respectively. An important observation is that the unitary $\mathcal{R}$ is
followed by the projective measurement in the computational basis
$\\{\ket{\mathbf{x}}\\}$, and thus the quantum gates that do not create
superposition before the final measurement can be equivalently moved to the
classical post-processing [48]. As such, the circuit in Fig. 3 (a) can be
effectively transformed to the one in Fig. 3 (b).

We can understand the structure of the circuit in Fig. 3 (b) as follows. Note
that the subspaces $\mathcal{V}$ and $\mathcal{V}^{\perp}$ depend on the
parity information of the input, and $CX_{1\rightarrow 2}$ followed by
measurement on the second qubit aims to obtain such parity information $c$.
This effectively projects (partially decoherence) the input into $\mathcal{V}$
or $\mathcal{V}^{\perp}$. Conditioned on the parity $c$, one adds a
_classical_ controlled-$H$ on the first qubit, thus creating superposition
only in $\mathcal{V}$ but not in $\mathcal{V}^{\perp}$. The measurement result
of the first qubit is $x_{1}$, and $x_{2}=x_{1}\oplus c$ by the classical
addition. One can extend the circuit in Fig. 3 (b) to the qudit case by
substituting the $CX$ and $H$ gates with their qudit analogs, and the detail
is left in Appendix C.

Figure 3: Circuit compilation of $\mathcal{R}$. (a) Compilation of
$\mathcal{R}$ in the single-qubit case. (b) Simplification of the circuit by
utilizing the classical post-processing. (c) The generalization to the
$n$-qubit case. Here, we take $n=4$, and assume that the parity check
$c^{1}=0$, and $c^{i}=1$ for $i\in\\{2,3,4\\}$. For the first qubit pair, the
quantum circuit is similar to the single-qubit case, but there is a joint
unitary consisting CNOTs among the last three qubits in the first replica. In
addition, as in the single-qubit case shown in (b), there are some classical
post-processings to calculate the final measurement label $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$.

The circuit compilation in the many-qubit case is non-trivial and indeed not a
direct extension from the single-qubit case, since
$\mathcal{R}\neq\bigotimes_{i}\mathcal{R}_{(i)}$, i.e., $\mathcal{R}$ cannot
be written in a tensor-product form. Here, $\mathcal{R}_{(i)}$ denotes the
local unitary on the $i$-th qubit pair. Note that the standard Bell
measurement nevertheless can be written into the qubit-wise tensor-product
form [15, 32], which highlights the significant difference between these two
bases. The main strategy of the following compilation is to add more
entangling operations inside the first and(or) second replicas, beyond the
product structure where the initial $R_{(i)}$ only supplies entanglement
across the two replicas.

Use $\mathbf{a}$ and $\mathbf{b}$ to label the $n$-qubit basis for the first
and second replicas, respectively. Without loss of generality, assume that the
first $n-l$ bits of $\mathbf{a}$ and $\mathbf{b}$ share the same string, and
other $l$ bits are different (opposite indeed). That is,
$\mathbf{a}=\mathbf{a}^{[n-l]}\mathbf{a}^{[l]}$ and
$\mathbf{b}=\mathbf{a}^{[n-l]}\overline{\mathbf{a}^{[l]}}$. Here, we use
$\mathbf{a}^{[n-l]}(\mathbf{a}^{[l]})$ to denote the restricted bit-string on
the subsystem, and $\overline{\mathbf{a}^{[l]}}$ denotes the bit-flip of
$\mathbf{a}^{[l]}$. By Eq. 8, one has

$\displaystyle\mathcal{R}\ket{\mathbf{a}^{[n-l]}\mathbf{a}^{[l]}}_{1}\ket{\mathbf{a}^{[n-l]}\overline{\mathbf{a}^{[l]}}}_{2}=\ket{\mathbf{a}^{[n-l]}}_{1}\ket{\mathbf{a}^{[n-l]}}_{2}\otimes\left\\{\begin{array}[]{lll}\\!\frac{1}{\sqrt{2}}\left(\ket{\overline{\mathbf{a}^{[l]}}}_{1}\ket{\mathbf{a}^{[l]}}_{2}+\ket{\mathbf{a}^{[l]}}_{1}\ket{\overline{\mathbf{a}^{[l]}}}_{2}\right)&\quad\mathbf{a}^{[l]}<\overline{\mathbf{a}^{[l]}},\\\\[6.02773pt]
\\!\frac{1}{\sqrt{2}}\left(\ket{\overline{\mathbf{a}^{[l]}}}_{1}\ket{\mathbf{a}^{[l]}}_{2}-\ket{\mathbf{a}^{[l]}}_{1}\ket{\overline{\mathbf{a}^{[l]}}}_{2}\right)&\quad\mathbf{a}^{[l]}>\overline{\mathbf{a}^{[l]}},\end{array}\right.$
(13)

where we just extract the common state on the first $n-l$ qubits on both
replicas; and the inequality $\mathbf{a}^{[l]}<\overline{\mathbf{a}^{[l]}}$
($\mathbf{a}^{[l]}>\overline{\mathbf{a}^{[l]}}$) means that the first bit of
$\mathbf{a}^{[l]}$ is 0 (1). It is not hard to see that the state on the last
$l$ qubits is a GHZ-like state, which inspire the following compilation.

Let us summarize the action of $\mathcal{R}$ in Eq. 13: if $a_{i}=b_{i}$,
$\mathcal{R}$ does nothing; and for the rest of qubits in $[l]$ with
$b_{i}=\overline{a_{i}}$, it creates a ‘large’ superposition in the subspace
$\mathcal{V}={\rm
span}\left(\ket{\mathbf{a}^{[l]}}\ket{\overline{\mathbf{a}^{[l]}}},\ket{\overline{\mathbf{a}^{[l]}}}\ket{\mathbf{a}^{[l]}}\right)$,
which degenerates to ${\rm span}\left(\ket{01},\ket{10}\right)$ in the single-
qubit case. By applying the technique of quantum error correction code [34],
we construct the quantum circuit as shown in Fig. 3 (c). Conditioned on the
parity check information on the second replica, we additionally add a joint
unitary on the qubits in the first replica where the parity $c^{i}=1$. Such
joint unitary is composed of CNOTs and a Hadamard gate, which are both
Clifford gates [34, 49] and can be realized in a fault-tolerant way in the
future [50]. We leave more details in Appendix C.

Some remarks on the structure of many-qubit quantum circuit in Fig. 3 (c) are
as follows. First, after the initial parity check which needs depth-1 circuit,
the remaining circuit depth equals the number of $c^{i}=1$, and thus the
maximal circuit depth could be $n+1$ in total. Second, one needs long-range
CNOTs which are realizable in systems like ion-trap [51] and optical-tweezer
platforms [52], but would be challenging for the ones with the nearest-
neighbour architecture, like superconducting-qubit system [20, 53]. The
compilation of the long-range gates in such systems would also increase the
circuit depth further. In the remaining sections, we develop advanced
measurement protocol beyond AFRS in Section .1 to significantly reduce the
circuit depth.

###  .4 Local-AFRS

In all the previous sections, we focus on constructing the estimator
$\widehat{\rho^{t}}$ with AFRS. Suppose the observable $O$ is local and acting
nontrivially on subsystem $A$, say $O=O_{A}\otimes\mathbb{I}_{\bar{A}}$. By
further developing AFRS for such local observables, we propose the local-AFRS
protocols, which can significantly reduce the difficulty of realizing the
entangling unitary $\mathcal{R}$ in Section .3, and in the same time keeps the
same estimation performance as AFRS. In particular, we show that one can
estimate $o_{t}=\tr(O\rho^{t})$, or more generally the reduced matrix
$\tr_{\bar{A}}(\rho^{t})$ on $A$ with greatly simplified quantum circuit, even
to constant depth when the qubit number in $A$ is $|A|=\mathcal{O}(1)$.

Figure 4: The auxiliary-free replica shadow protocol of estimating
$\tr(O\rho^{t})$ for some local observables. (a) In local-AFRS, the random
unitary evolution takes the product form $V=V_{A}\otimes V_{\bar{A}}$ across
the subsystems $A$ and $\bar{A}$ for all $t$ replicas. The joint unitary
before the final computational basis measurement can also take the product
form, $\mathcal{R}_{A}\otimes\mathcal{R}_{\bar{A}}$, where
$\mathcal{R}_{A(\bar{A})}$ is the entangling unitary restricted on $A$
$(\bar{A})$ for all replicas. The circuit depth of $\mathcal{R}_{A(\bar{A})}$
is $|A|+1$ $(|\bar{A}|+1)$ for $t=2$. With the randomized measurement results,
one can follow Eq. 17 in Theorem 2 to build the unbiased estimator of
$\rho^{t}$ reduced on subsystem $A(\bar{A})$ as
$\widehat{\tr_{\bar{A}}(\rho^{t})}$
$\Big{(}\widehat{\tr_{A}(\rho^{t})}\Big{)}$, which can be used to estimate any
$O_{A}$ and $O_{\bar{A}}$. (b) If one only cares the observables on $A$, the
random unitary on $\bar{A}$ can take identity, i.e.,
$V_{\bar{A}}=\mathbb{I}_{{\bar{A}}}$. And $\mathcal{R}_{\bar{A}}$ can be
further substituted by the qubit-wise circuit
$\bigotimes_{i\in\bar{A}}\mathcal{R}_{(i)}$. Here, $\mathcal{R}_{(i)}$ is the
entangling unitary in single-qubit case among $t$ replicas. The depth of the
circuit on $\bar{A}$ is at most $2$ for $t=2$.

By respecting the locality property, we let the random unitary in the product
form $V=V_{A}\otimes V_{\bar{A}}$. Similar as the main framework of AFRS in
Section .1, here one needs to effectively generate the _fake_ but now
_marginal_ conditional probability
$\Pr(\mathbf{b}^{A}|V_{A})=\bra{\mathbf{b}^{A}}V_{A}(\tr_{\bar{A}}(\rho^{t}))\ket{\mathbf{b}^{A}}$,
where $\mathbf{b}^{A}$ labels the bit-string restricted on $A$. It is not hard
to check that 1 in 1 can be applied to the local case here, but with the
observable $Q_{\mathbf{b}}$ in 1 updated to

$\displaystyle Q^{\mathrm{local}}(\mathbf{b}^{A})$
$\displaystyle=\left[t^{-1}\sum_{i}\ket{\mathbf{b}^{A}}\bra{\mathbf{b}^{A}}_{i}\otimes\mathbb{I}^{\otimes(t-1)}_{A}\right]\otimes\mathbb{I}^{\otimes
t}_{\bar{A}}$ (14) $\displaystyle=Q_{\mathbf{b}^{A}}\otimes\mathbb{I}^{\otimes
t}_{\bar{A}}.$

Here, $Q_{\mathbf{b}^{A}}$ is symmetric on all $t$ replicas, but restricted on
$A$; the complement subsystem $\bar{A}$ takes the identity
$\mathbb{I}_{\bar{A}}$ for all $t$ replicas. We explain in detail and prove
such observation in Section D.1.

Remember that the key step to develop AFRS in Section .1 is to find the common
basis of both $S_{t}$ and $Q_{\mathbf{b}}$. Note that now the updated
$Q^{\mathrm{local}}(\mathbf{b}^{A})$ in Eq. 14 is in the tensor-product form.
Importantly, $S_{t}$ can always be factorized to, $S_{t}=S_{t}^{A}\otimes
S_{t}^{\bar{A}}$, and even further to the qubit level on $t$ replicas. In this
way, the basis
$\mathcal{R}^{{\dagger}}_{A}\otimes\mathcal{R}^{{\dagger}}_{\bar{A}}\ket{\mathbf{x}^{A},\mathbf{x}^{\bar{A}}}=\ket{\psi_{\mathbf{x}^{A}}}\ket{\psi_{\mathbf{x}^{\bar{A}}}}$
are the common basis for both $Q^{\mathrm{local}}(\mathbf{b}^{A})$ and
$S_{t}$, with
$\mathbf{x}^{A}=\\{\mathbf{x}_{1}^{A}\mathbf{x}_{2}^{A}\cdots\mathbf{x}_{t}^{A}\\}$
for the bit-strings restricted on $A$ from all replicas. Thus one can use such
basis in the measurement, and the corresponding circuit is shown in Fig. 4
(a).

Accordingly, for local-AFRS here, the function $f(\mathbf{x})$ and the mapping
$\mathcal{G}_{\mathbf{x}\to\mathbf{b}}$ in 1 are updated respectively to

$\displaystyle f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})$
$\displaystyle=\bra{\psi_{\mathbf{x}^{A}}}\bra{\psi_{\mathbf{x}^{\bar{A}}}}S_{t}^{A}\otimes
S_{t}^{\bar{A}}\ket{\psi_{\mathbf{x}^{A}}}\ket{\psi_{\mathbf{x}^{\bar{A}}}}$
$\displaystyle=f(\mathbf{x}^{A})f(\mathbf{x}^{\bar{A}});$ (15)
$\displaystyle\mathcal{G}_{\mathbf{x}^{A}\to\mathbf{b}^{A}}$
$\displaystyle:\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})=\bra{\psi_{\mathbf{x}^{A}}}Q_{\mathbf{b}^{A}}\ket{\psi_{\mathbf{x}^{A}}},$
(16)

with the mapping in Eq. 16 restricted on $A$. With such construction, we have
the estimator and its performance of local-AFRS as follows.

######  Theorem 2.

Suppose one conducts the circuit in Fig. 4 (a) for once (single-shot, $M=1$),
the unbiased estimator of $\rho^{t}$ reduced on subsystem $A$ shows

$\displaystyle\widehat{\tr_{\bar{A}}(\rho^{t})}=\mathrm{Re}(f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}}))\mathcal{M}_{A}^{-1}\left(V_{A}^{{\dagger}}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right).$
(17)

Here, the inverse channel $\mathcal{M}_{A}^{-1}$ on $A$ depends on the random
unitary ensemble $V_{A}\in\mathcal{E}_{A}$;
$\mathbf{x}=\\{\mathbf{x}^{A},\mathbf{x}^{\bar{A}}\\}$ represents the actual
measurement outcome from the quantum circuit in Fig. 4 (a); the function $f$
is given in Eq. 15; and $\mathbf{b}^{A}$ is the mapping result by Eq. 16.

The local observable $O=O_{A}\otimes\mathbb{I}_{\bar{A}}$ on $\rho^{t}$ can be
estimated by
$\widehat{o_{t}}=\tr\\!\left(\widehat{\tr_{\bar{A}}(\rho^{t})}\,O_{A}\right)$,
and its variance shows

$\displaystyle\mathrm{Var}(\widehat{o_{t}})\leq\left\|O_{A,0}\right\|_{\mathrm{sh},\mathcal{E}_{A}}^{2}+\|O\|_{\infty}^{2},$
(18)

where $O_{A,0}$ is the traceless part of $O_{A}$; and
$\|\cdot\|_{\mathrm{sh},\mathcal{E}_{A}}$ is the shadow norm depending on the
unitary ensemble $\mathcal{E}_{A}$.

Theorem 1 can be seen as a special case of Theorem 2, with the subsystem being
the whole system, say $A=[n]$, and the proof of Theorem 2 is left in Section
D.2.

If focusing on the estimation of $O_{A}$, one only needs the unitary
$V=V_{A}\otimes\mathbb{I}_{\bar{A}}$ in the random evolution, that is, the
unitary on $\bar{A}$ can take the identity. In addition, the entangling
operation
$\mathcal{R}_{A}\otimes\mathcal{R}_{\bar{A}}\rightarrow\mathcal{R}_{A}\otimes\bigotimes_{i\in\bar{A}}\mathcal{R}_{(i)}$
can be further factorized to qubit-pair level on $\bar{A}$, as shown in Fig. 4
(b). For the important case $t=2$, the circuit depth is at most $|A|+1$
(mainly used to realize $\mathcal{R}_{A}$), being usually some constant number
if $|A|=\mathcal{O}(1)$. Consequently, compared with AFRS, local-AFRS
significantly reduces the circuit depth from $n+1$ to $|A|+1$ for $t=2$, but
keeps the estimation performance almost the same as AFRS, indicated by Eq. 18.
Moreover, suppose one considers the extreme case $O=\mathbb{I}$ and
$A=\emptyset$. Then the estimator in Eq. 17 becomes that for the $t$-th moment
showing,

$\displaystyle\widehat{P_{t}}:=\widehat{\tr(\rho^{t})}=\mathrm{Re}\left(\prod_{i=1}^{n}f(\mathbf{x}^{i})\right),$
(19)

where $\mathbf{x}^{i}$ is the measurement result for the $i$-th qubit on all
$t$ replicas after the local rotation $\bigotimes_{i=1}^{n}\mathcal{R}_{(i)}$.
By Eq. 18 and its proof in Section D.2, the variance shows
$\mathrm{Var}(\widehat{P_{t}})\leq 1-\left[\tr(\rho^{t})\right]^{2}$. We give
the corresponding proof and more explanations of such moment function in
Section D.3.

These facts make local-AFRS more feasible in practical quantum platforms,
which are further manifested numerically in Fig. 2 (a). By using OS, AFRS and
local-AFRS to estimate $\tr(O\rho^{2})$ with $O$ being some local Pauli
operator, we show that the performance of local-AFRS is almost the same as
AFRS, with an exponential advantage compared to OS with respective to the
estimation-error scaling to the system size [19, 33].

One can pursue such circuit complexity reduction further by considering the
problem of estimating a collection of local observables
$\mathcal{W}=\\{O_{l}\\}_{l=1}^{L}$, which may have different local supports.
Divide the observable set as $\mathcal{W}=\cup_{k=1}^{K}\mathcal{W}_{k}$ with
each of the subset satisfying

$\displaystyle\mathcal{W}_{k}=\\{O_{l}|\ \forall\\{l,l^{\prime}\\},$
$\displaystyle\operatorname{\mathrm{supp}}(O_{l})\cap\operatorname{\mathrm{supp}}(O_{l^{\prime}})=\emptyset$
(20)
$\displaystyle\mathrm{or}\operatorname{\mathrm{supp}}(O_{l})\subseteq\operatorname{\mathrm{supp}}(O_{l^{\prime}})\\}.$

For each observable subset $\mathcal{W}_{k}$, one can apply the entangling
unitary $U_{k}=\bigotimes_{j}\mathcal{R}_{A_{j}}$ before the final projective
measurement, such that $\\{A_{j}\\}$ can cover the support of any
$O_{l}\in\mathcal{W}_{k}$. In this way, one needs to apply $K$ different
$U_{k}$ in the experimental stage one at a time. It is clear the depth of the
circuit is at most $\max_{l}\operatorname{\mathrm{supp}}(O_{l})+1$, and we
show the performance for such strategy as follows, which inherits the nature
of shadow estimation and maintains the sampling advantage of AFRS.

######  Corollary 1.

Given a collection of local observables $\mathcal{W}=\\{O_{l}\\}_{l=1}^{L}$,
divide them into $K$ subsets $\\{\mathcal{W}_{k}\\}_{k=1}^{K}$ satisfying Eq.
20. Suppose one applies local-AFRS to estimate them. Then the repeating
measurement time or the size of the shadow set,

$\displaystyle
M=\frac{68K}{\epsilon^{2}}\max_{O_{l}\in\mathcal{W}}\mathrm{Var}\Big{(}\widehat{o^{(l)}_{t}}\Big{)}\log\left(\frac{2L}{\delta}\right),$
(21)

is sufficient to guarantee that the estimation error of each
$O_{l}\in\mathcal{W}$ is less than $\epsilon$, with probability at least
$1-\delta$. Note that $\widehat{o^{(l)}_{t}}$ is the single-shot estimator of
the observable $O_{l}$ on $\rho^{t}$.

Here, we use the median-of-means post-processing method to the shadow sets
[19]; see Appendix F for the proof and more details. Suppose all
$\operatorname{\mathrm{supp}}(O_{l})=\mathcal{O}(1)$, the variance in Eq. 21
for each observable can thus be bounded by some constant using Eq. 18 in
Theorem 2. In this way, under a constant estimation error and confidence
level, the measurement time in Eq. 21 scales like $M=\mathcal{O}(K\log(L))$.
The linear scaling with $K$ is due to the implementation of $K$ different
entangling unitary in the estimation, and the logarithmic scaling with $L$ is
on account of the shadow nature of local-AFRS, which multiplexes the shadow
snapshots in the post-processing. For instance, suppose $K=\mathcal{O}(1)$ and
$L=\mathcal{O}(n)$, one has $M=\mathcal{O}(\log(n))$ compared to previous
methods with $M=\mathcal{O}(n)$ [15, 16]. We clarify such advantages of local-
AFRS more in-depth with the virtual distillation task in the following
section.

Figure 5: Estimation performance of OS, AFRS and local-AFRS protocols. The
quantum state is a noisy $5$-qubit GHZ state
$\rho=0.7\ket{\text{GHZ}}\bra{\text{GHZ}}+0.3{\mathbb{I}}_{d}/d$. Local
Clifford ensemble with $V\in\mathcal{E}_{\mathrm{LCl}}$ is utilized for
evaluating $o_{2}=\tr(O\rho^{2})$ in (a), and $\langle
O\rangle_{\mathrm{VD}}^{(t=2)}=\tr(O\rho^{2})/\tr(\rho^{2})$ in (b). In (a),
we set $M=500$ for AFRS and local-AFRS protocols, and $M=1000$ for OS
protocol, thus the comparison is carried out under the same sample number of
$\rho$ among different protocols. In (b), the local observable is set as
$O=Z_{1}Z_{2}$, and the black dashed line represents the theoretical result of
$\langle O\rangle_{\mathrm{VD}}^{(t=2)}=0.994$.

###  .5 Application in virtual distillation

Virtual distillation (VD) aims to purify the mixed state $\rho$ using its
$t$-th order ‘virtual’ state $\rho^{(t)}_{\text{VD}}=\rho^{t}/\tr(\rho^{t})$
[14, 15, 16]. In general, the transformation
$\rho\rightarrow\rho^{(t)}_{\text{VD}}$ is not a physical operation, and one
can alternatively make up the expectation value of some observable $O$, say
$\langle
O\rangle^{(t)}_{\text{VD}}=\tr(O\rho^{(t)}_{\text{VD}})=\tr(O\rho^{t})/\tr(\rho^{t})$,
by some quantum measurements and post-processings. VD can be applied to
algorithmic cooling of many-body systems [14], and to quantum error mitigation
protocols by amplifying the dominant state vector [15, 16, 31].

Here, we demonstrate the strengths of local-AFRS described in Section .4 in
such VD task. Consider the Hamiltonian of the transverse-field Ising model
$H=-J\sum_{j}Z_{j}Z_{j+1}-h\sum_{j}X_{j}$, and we would like to estimate
$\langle O\rangle^{(t)}_{\text{VD}}$ for all Pauli operators in the summation,
say $O\in\\{Z_{j}Z_{j+1},X_{j}\\}$. Note that the moment function in the
denominator $\tr(\rho^{t})$ can be estimated by just taking $O=\mathbb{I}$,
and thus the essential part is the numerator $o_{t}=\tr(O\rho^{t})$ for
various $O$. By Corollary 1, we only need to choose $K=2$ different entangling
unitary to accomplish the estimation. That is, after the local random Clifford
evolution $V=\bigotimes_{i=1}^{n}V_{i}$ on each replica, one applies the
following two entangling unitaries one at a time (with $n$ being odd for
simplicity),

$\displaystyle U_{\text{odd}}=\bigotimes_{j=1}^{\lfloor
n/2\rfloor}\mathcal{R}_{2j-1,2j}\otimes\mathcal{R}_{n},$ (22) $\displaystyle
U_{\text{even}}=\mathcal{R}_{1}\otimes\bigotimes_{j=1}^{{\lfloor
n/2\rfloor}}\mathcal{R}_{2j,2j+1},$

to measure all $t$ replicas of $\rho$. For $\\{Z_{j}Z_{j+1}\\}$ with an
odd/even $j$, one post-processes the measurement results from
$U_{\text{odd}/\text{even}}$, respectively. For $\\{X_{i}\\}$ and
$\mathbb{I}$, both measurement patterns could be used, as they cover the
supports of the observables. In this way, the sampling time is just
$M=\mathcal{O}(\log(n))$ by Corollary 1. Note that for the case $t=2$, the
total circuit depth is at most $3$.

Some remarks on the advantages with respect to the previous methods are
illustrated as follows. First, compared with Ref. [16], local-AFRS does not
need any auxiliary qubit. Most importantly, compared to a constant-depth
circuit here, the generalized controlled-swap operation among all replicas
used there would make the depth of the circuit being $\mathcal{O}(n)$. This
fact indicates that our protocol is more feasible on current and near-future
platforms. Second, the auxiliary-free protocol based on Bell-like measurements
given in Ref. [15], mainly works for single-qubit observables, and even the
generalization to two-qubit observable needs to find some specific rotation
unitary case-by-case before the final measurement [31]. Thus, it is quite hard
and not clear to generalize to $k$-local ones. In addition, the sampling time
there scales as $\mathcal{O}(n)$, which is exponentially worse than
$\mathcal{O}(\log(n))$ presented here. Last, compared with OS [19, 33], our
protocol harnesses the power of joint unitary and thus shows an exponential
advantage on sampling time as already shown in Section .4.

In Fig. 5, we also numerically manifest the advantages of local-AFRS for
estimating Pauli observables $O$ and in the context of VD, compared to the OS
protocol. Fig. 5 (a) shows the estimation error of $\tr(O\rho^{2})$ for OS,
AFRS and local-AFRS with various observables
$O=\\{Z_{j}Z_{j+1},X_{j},{\mathbb{I}}\\}$ in the Ising Hamiltonian. The
performances of AFRS and local-AFRS protocols are similar, and much better
than that of OS. Their performances depend on the size of the support of the
observables, which are consistent with the variance results predicted in
Theorem 1 and Theorem 2. Fig. 5 (b) further compares the VD results of
$\langle O\rangle^{(t=2)}_{\text{VD}}$ for OS and local-AFRS under the same
sampling number of the quantum state $\rho$. It is clear that local-AFRS
converges faster, and the VD result is much closer to the theoretical one than
OS, which demonstrates the strength of local-AFRS in VD.

## Discussions

In this study, we proposed AFRS framework for efficiently estimating general
nonlinear functions of quantum states, by utilizing the power of entangling
measurements on a few state replicas without using any auxiliary qubits.
Different from traditional destructive swap test, the designed entangling
measurement together with the foregoing random unitary evolution and classical
post-mapping and processing, can surprisingly return an unbiased estimator of
complex nonlinear functional. Inheriting the mindset of shadow estimation, the
measurement data can be further reused for various observables and functions.
Moreover, with the improved local-AFRS protocol, one can realize the
estimation of local observables with even constant-depth quantum circuits,
which facilitates its large-scale realization even on the current quantum
platforms. AFRS would definitely serve as a significant component, similar as
the famous swap test, in quantum information processing, which would inspire
many developments and applications from quantum algorithm design to quantum
many-body physics.

The questions and solutions endorsed in this work should be a beginning rather
than the conclusion of the study. First, it is intriguing to extend the
current framework to fermion and boson systems in quantum simulation [54, 55,
56], and to quantum channel and dynamics [57, 58, 59] as well, where the key
ingredient should be the fine design of entangling measurements in all
scenarios. Second, for the many-qubit scenario considered in the current work,
it is interesting and important to find the corresponding circuit compilation
when the replica number $t\geq 3$ even for small-size subsystems, which would
be helpful to advance the estimation of higher-order functions [60, 61, 62].
Third, AFRS can be further integrated with other shadow schemes, such as
hybrid shadow [21], conjugate shadow [63], multi-shot shadow [46, 47] (some
preliminary results are shown in Appendix E), and shallow shadow [64, 65, 66],
to enhance its abilities further in the estimation. In particular, AFRS would
play a significant role in the principal component analysis [67, 68, 69].
Fourth, considering the quantum error mitigation [70, 71] and the noise
aspect, starting from the VD demonstrated here, we expect AFRS together with
other mitigation techniques and protocols would advance the development of
quantum error mitigation [72, 73] and the robustness of the measurement scheme
itself [74, 75, 76, 77]. Finally, as mentioned before, a number of important
quantities like quantum Fisher information [43, 44] and out-of-time-ordered
correlators [11, 13], can be expressed as or approximated by nonlinear
functions considered in this work, we expect the realization of AFRS to
measure them on large-scale quantum platforms for practical applications.

## Acknowledgements

Q.L.$\&$Y.Z. acknowledge the support from the National Natural Science
Foundation of China (NSFC) (Grant No. 12205048), Innovation Program for
Quantum Science and Technology (Grant No. 2021ZD0302000), and the start-up
funding of Fudan University. Z.L.$\&$H.Z. acknowledge the support from NSFC
(Grant No. 92165109), National Key Research and Development Program of China
(Grant No. 2022YFA1404204), and Shanghai Municipal Science and Technology
Major Project (Grant No. 2019SHZDZX01). X.Y. acknowledges the support from
NSFC (Grant No. 12175003 and No. 12361161602), NSAF (Grant No. U2330201), the
Innovation Program for Quantum Science and Technology (Grant No.
2023ZD0300200).

## References

  * Ekert  _et al._ [2002] A. K. Ekert, C. M. Alves, D. K. L. Oi, M. Horodecki, P. Horodecki, and L. C. Kwek, Direct estimations of linear and nonlinear functionals of a quantum state, Phys. Rev. Lett. 88, 217901 (2002).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * Johri  _et al._ [2017] S. Johri, D. S. Steiger, and M. Troyer, Entanglement spectroscopy on a quantum computer, Phys. Rev. B 96, 195136 (2017).
  * Bruni [2004] T. A. Bruni, Measurimg polynomial functions of states, Quantum Info. Comput. 4, 401–408 (2004).
  * Pichler  _et al._ [2016] H. Pichler, G. Zhu, A. Seif, P. Zoller, and M. Hafezi, Measurement protocol for the entanglement spectrum of cold atoms, Phys. Rev. X 6, 041033 (2016).
  * Elben  _et al._ [2020a] A. Elben, R. Kueng, H.-Y. R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller, and B. Vermersch, Mixed-state entanglement from local randomized measurements, Phys. Rev. Lett. 125, 200501 (2020a).
  * Zhou  _et al._ [2020] Y. Zhou, P. Zeng, and Z. Liu, Single-copies estimation of entanglement negativity, Phys. Rev. Lett. 125, 200502 (2020).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. De Wolf, Quantum fingerprinting, Phys. Rev. Lett. 87, 167902 (2001).
  * Horodecki and Ekert [2002] P. Horodecki and A. Ekert, Method for direct detection of quantum entanglement, Phys. Rev. Lett. 89, 127902 (2002).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).
  * Vermersch  _et al._ [2019] B. Vermersch, A. Elben, L. M. Sieberer, N. Y. Yao, and P. Zoller, Probing scrambling using statistical correlations between randomized measurements, Phys. Rev. X 9, 021061 (2019).
  * Elben  _et al._ [2020b] A. Elben, J. Yu, G. Zhu, M. Hafezi, F. Pollmann, P. Zoller, and B. Vermersch, Many-body topological invariants from randomized measurements in synthetic quantum matter, Sci. Adv. 6, eaaz3666 (2020b).
  * Garcia  _et al._ [2021] R. J. Garcia, Y. Zhou, and A. Jaffe, Quantum scrambling with classical shadows, Phys. Rev. Res. 3, 033155 (2021).
  * Cotler  _et al._ [2019] J. Cotler, S. Choi, A. Lukin, H. Gharibyan, T. Grover, M. E. Tai, M. Rispoli, R. Schittko, P. M. Preiss, A. M. Kaufman, M. Greiner, H. Pichler, and P. Hayden, Quantum virtual cooling, Phys. Rev. X 9, 031013 (2019).
  * Huggins  _et al._ [2021] W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, Virtual distillation for quantum error mitigation, Phys. Rev. X 11, 041036 (2021).
  * Koczor [2021] B. Koczor, Exponential error suppression for near-term quantum devices, Phys. Rev. X 11, 031057 (2021).
  * Elben  _et al._ [2023] A. Elben, S. T. Flammia, H.-Y. Huang, R. Kueng, J. Preskill, B. Vermersch, and P. Zoller, The randomized measurement toolbox, Nat. Rev. Phys. 5, 9 (2023).
  * Aaronson [2019] S. Aaronson, Shadow tomography of quantum states, SIAM Journal on Computing 49, STOC18 (2019).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nat. Phys. 16, 1050 (2020).
  * Huang  _et al._ [2022] H.-Y. Huang, M. Broughton, J. Cotler, S. Chen, J. Li, M. Mohseni, H. Neven, R. Babbush, R. Kueng, J. Preskill, _et al._ , Quantum advantage in learning from experiments, Science 376, 1182 (2022).
  * Zhou and Liu [2024] Y. Zhou and Z. Liu, A hybrid framework for estimating nonlinear functions of quantum states, npj Quantum Inf. 10, 62 (2024).
  * Quek  _et al._ [2024] Y. Quek, E. Kaur, and M. M. Wilde, Multivariate trace estimation in constant quantum depth, Quantum 8, 1220 (2024).
  * Patel  _et al._ [2016] R. B. Patel, J. Ho, F. Ferreyrol, T. C. Ralph, and G. J. Pryde, A quantum Fredkin gate, Sci. Adv. 2, e1501531 (2016).
  * Gao  _et al._ [2019] Y. Y. Gao, B. J. Lester, K. S. Chou, L. Frunzio, M. H. Devoret, L. Jiang, S. Girvin, and R. J. Schoelkopf, Entanglement of bosonic modes through an engineered exchange interaction, Nature 566, 509 (2019).
  * Zhang  _et al._ [2019] K. Zhang, J. Thompson, X. Zhang, Y. Shen, Y. Lu, S. Zhang, J. Ma, V. Vedral, M. Gu, and K. Kim, Modular quantum computation in a trapped ion system, Nat. Commun. 10, 4692 (2019).
  * Gan  _et al._ [2020] H. Gan, G. Maslennikov, K.-W. Tseng, C. Nguyen, and D. Matsukevich, Hybrid quantum computing with conditional beam splitter gate in trapped ion system, Phys. Rev. Lett. 124, 170502 (2020).
  * Peng  _et al._ [2024] X.-J. Peng, Q. Liu, L. Liu, T. Zhang, Y. Zhou, and H. Lu, Experimental hybrid shadow tomography and distillation, arXiv preprint arXiv:2404.11850 (2024).
  * Subaşı  _et al._ [2019] Y. Subaşı, L. Cincio, and P. J. Coles, Entanglement spectroscopy with a depth-two quantum circuit, J. Phys. A: Math. Theor. 52, 044001 (2019).
  * Yirka and Subaşı [2021] J. Yirka and Y. Subaşı, Qubit-efficient entanglement spectroscopy using qubit resets, Quantum 5, 535 (2021).
  * Bluvstein  _et al._ [2024] D. Bluvstein, S. J. Evered, A. A. Geim, S. H. Li, H. Zhou, T. Manovitz, S. Ebadi, M. Cain, M. Kalinowski, D. Hangleiter, _et al._ , Logical quantum processor based on reconfigurable atom arrays, Nature 626, 58 (2024).
  * O’Brien  _et al._ [2023] T. E. O’Brien, G. Anselmetti, F. Gkritsis, V. Elfving, S. Polla, W. J. Huggins, O. Oumarou, K. Kechedzhi, D. Abanin, R. Acharya, _et al._ , Purification-based quantum error mitigation of pair-correlated electron simulations, Nat. Phys. 19, 1787 (2023).
  * Subaşı [2022] Y. Subaşı, Ancilla-free continuous-variable SWAP test, Quantum 6, 800 (2022).
  * Seif  _et al._ [2023] A. Seif, Z.-P. Cian, S. Zhou, S. Chen, and L. Jiang, Shadow distillation: Quantum error mitigation with classical shadows for near-term quantum processors, PRX Quantum 4, 010303 (2023).
  * Gottesman [2002] D. Gottesman, An introduction to quantum error correction, in  _Proceedings of Symposia in Applied Mathematics_, Vol. 58 (2002) pp. 221–236.
  * Hu and You [2022] H.-Y. Hu and Y.-Z. You, Hamiltonian-driven shadow tomography of quantum states, Phys. Rev. Res. 4, 013054 (2022).
  * Hu  _et al._ [2023a] H.-Y. Hu, S. Choi, and Y.-Z. You, Classical shadow tomography with locally scrambled quantum dynamics, Phys. Rev. Res. 5, 023027 (2023a).
  * Ohliger  _et al._ [2013] M. Ohliger, V. Nesme, and J. Eisert, Efficient and feasible state tomography of quantum many-body systems, New J. Phys. 15, 015024 (2013).
  * Chen  _et al._ [2022] S. Chen, J. Cotler, H.-Y. Huang, and J. Li, Exponential separations between learning with and without quantum memory, in  _2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)_ (IEEE, 2022) pp. 574–585.
  * Elben  _et al._ [2020c] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, Cross-platform verification of intermediate scale quantum devices, Phys. Rev. Lett. 124, 010504 (2020c).
  * Liu  _et al._ [2022a] Z. Liu, P. Zeng, Y. Zhou, and M. Gu, Characterizing correlation within multipartite quantum systems via local randomized measurements, Phys. Rev. A 105, 022407 (2022a).
  * van Enk and Beenakker [2012] S. J. van Enk and C. W. J. Beenakker, Measuring $\mathrm{Tr}{\rho}^{n}$ on single copies of $\rho$ using random measurements, Phys. Rev. Lett. 108, 110503 (2012).
  * Brydges  _et al._ [2019] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Rényi entanglement entropy via randomized measurements, Science 364, 260 (2019).
  * Rath  _et al._ [2021] A. Rath, C. Branciard, A. Minguzzi, and B. Vermersch, Quantum Fisher information from randomized measurements, Phys. Rev. Lett. 127, 260501 (2021).
  * Yu  _et al._ [2021a] M. Yu, D. Li, J. Wang, Y. Chu, P. Yang, M. Gong, N. Goldman, and J. Cai, Experimental estimation of the quantum Fisher information from randomized measurements, Phys. Rev. Res. 3, 043122 (2021a).
  * Cian  _et al._ [2021] Z.-P. Cian, H. Dehghani, A. Elben, B. Vermersch, G. Zhu, M. Barkeshli, P. Zoller, and M. Hafezi, Many-body chern number from statistical correlations of randomized measurements, Phys. Rev. Lett. 126, 050501 (2021).
  * Helsen and Walter [2023] J. Helsen and M. Walter, Thrifty shadow estimation: Reusing quantum circuits and bounding tails, Phys. Rev. Lett. 131, 240602 (2023).
  * Zhou and Liu [2023] Y. Zhou and Q. Liu, Performance analysis of multi-shot shadow estimation, Quantum 7, 1044 (2023).
  * Liu  _et al._ [2022b] G. Liu, X. Zhang, and X. Ma, Classically replaceable operations, Quantum 6, 845 (2022b).
  * Ozols [2008] M. Ozols, Clifford group, Essays at University of Waterloo, Spring (2008).
  * Campbell  _et al._ [2017] E. T. Campbell, B. M. Terhal, and C. Vuillot, Roads towards fault-tolerant universal quantum computation, Nature 549, 172 (2017).
  * Figgatt  _et al._ [2019] C. Figgatt, A. Ostrander, N. M. Linke, K. A. Landsman, D. Zhu, D. Maslov, and C. Monroe, Parallel entangling operations on a universal ion-trap quantum computer, Nature 572, 368 (2019).
  * Evered  _et al._ [2023] S. J. Evered, D. Bluvstein, M. Kalinowski, S. Ebadi, T. Manovitz, H. Zhou, S. H. Li, A. A. Geim, T. T. Wang, N. Maskara, _et al._ , High-fidelity parallel entangling gates on a neutral atom quantum computer, Nature 622, 268–272 (2023).
  * Cao  _et al._ [2023] S. Cao, B. Wu, F. Chen, M. Gong, Y. Wu, Y. Ye, C. Zha, H. Qian, C. Ying, S. Guo, _et al._ , Generation of genuine entanglement up to 51 superconducting qubits, Nature 619, 738 (2023).
  * Gandhari  _et al._ [2022] S. Gandhari, V. V. Albert, T. Gerrits, J. M. Taylor, and M. J. Gullans, Continuous-variable shadow tomography, arXiv preprint arXiv:2211.05149 10.48550/arXiv.2211.05149 (2022).
  * Becker  _et al._ [2024] S. Becker, N. Datta, L. Lami, and C. Rouzé, Classical shadow tomography for continuous variables quantum systems, IEEE Trans. Inf. Theory (2024).
  * Zhao  _et al._ [2021] A. Zhao, N. C. Rubin, and A. Miyake, Fermionic partial tomography via classical shadows, Phys. Rev. Lett. 127, 110504 (2021).
  * Helsen  _et al._ [2023] J. Helsen, M. Ioannou, J. Kitzinger, E. Onorati, A. Werner, J. Eisert, and I. Roth, Shadow estimation of gate-set properties from random sequences, Nat. Commun. 14, 5039 (2023).
  * Kunjummen  _et al._ [2023] J. Kunjummen, M. C. Tran, D. Carney, and J. M. Taylor, Shadow process tomography of quantum channels, Phys. Rev. A 107, 042403 (2023).
  * Liu  _et al._ [2024] Z. Liu, X. Zhang, Y.-Y. Fei, and Z. Cai, Virtual channel purification, arXiv preprint arXiv:2402.07866 (2024).
  * Ketterer  _et al._ [2019] A. Ketterer, N. Wyderka, and O. Gühne, Characterizing multipartite entanglement with moments of random correlations, Phys. Rev. Lett. 122, 120505 (2019).
  * Yu  _et al._ [2021b] X.-D. Yu, S. Imai, and O. Gühne, Optimal entanglement certification from moments of the partial transpose, Phys. Rev. Lett. 127, 060504 (2021b).
  * Liu  _et al._ [2022c] Z. Liu, Y. Tang, H. Dai, P. Liu, S. Chen, and X. Ma, Detecting entanglement in quantum many-body systems via permutation moments, Phys. Rev. Lett. 129, 260501 (2022c).
  * King  _et al._ [2024] R. King, K. Wan, and J. McClean, Exponential learning advantages with conjugate states and minimal quantum memory, arXiv preprint arXiv:2403.03469 (2024).
  * Hu  _et al._ [2023b] H.-Y. Hu, S. Choi, and Y.-Z. You, Classical shadow tomography with locally scrambled quantum dynamics, Phys. Rev. Res. 5, 023027 (2023b).
  * Bertoni  _et al._ [2024] C. Bertoni, J. Haferkamp, M. Hinsche, M. Ioannou, J. Eisert, and H. Pashayan, Shallow shadows: Expectation estimation using low-depth random Clifford circuits, Phys. Rev. Lett. 133, 020602 (2024).
  * Ippoliti  _et al._ [2023] M. Ippoliti, Y. Li, T. Rakovszky, and V. Khemani, Operator relaxation and the optimal depth of classical shadows, Phys. Rev. Lett. 130, 230403 (2023).
  * Lloyd  _et al._ [2014] S. Lloyd, M. Mohseni, and P. Rebentrost, Quantum principal component analysis, Nat. Phys. 10, 631 (2014).
  * Hakoshima  _et al._ [2023] H. Hakoshima, S. Endo, K. Yamamoto, Y. Matsuzaki, and N. Yoshioka, Localized virtual purification, arXiv preprint arXiv:2308.13500 (2023).
  * Grier  _et al._ [2024] D. Grier, H. Pashayan, and L. Schaeffer, Principal eigenstate classical shadows, arXiv preprint arXiv:2405.13939 (2024).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, Journal of the Physical Society of Japan 90, 032001 (2021).
  * Cai  _et al._ [2023] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum error mitigation, Rev. Mod. Phys. 95, 045005 (2023).
  * Huo and Li [2022] M. Huo and Y. Li, Dual-state purification for practical quantum error mitigation, Phys. Rev. A 105, 022427 (2022).
  * Cai [2021] Z. Cai, Quantum error mitigation using symmetry expansion, Quantum 5, 548 (2021).
  * Chen  _et al._ [2021] S. Chen, W. Yu, P. Zeng, and S. T. Flammia, Robust shadow estimation, PRX Quantum 2, 030348 (2021).
  * Koh and Grewal [2022] D. E. Koh and S. Grewal, Classical shadows with noise, Quantum 6, 776 (2022).
  * Wu and Koh [2024] B. Wu and D. E. Koh, Error-mitigated fermionic classical shadows on noisy quantum devices, npj Quantum Inf. 10, 39 (2024).
  * Jnane  _et al._ [2024] H. Jnane, J. Steinberg, Z. Cai, H. C. Nguyen, and B. Koczor, Quantum error mitigated classical shadows, PRX Quantum 5, 010324 (2024).
  * Jerrum  _et al._ [1986] M. R. Jerrum, L. G. Valiant, and V. V. Vazirani, Random generation of combinatorial structures from a uniform distribution, Theor. Comput. Sci. 43, 169 (1986).

###### Contents

  1. .1 Framework of AFRS
  2. .2 The matrix form of $\mathcal{R}$
  3. .3 Circuit compilation of $\mathcal{R}$
  4. .4 Local-AFRS
  5. .5 Application in virtual distillation
  6. A Proofs of 1 and Theorem 1
    1. A.1 Proof of 1
    2. A.2 Proof of Theorem 1
      1. A.2.1 Unbiasedness of the estimator $\widehat{\rho^{t}}$
      2. A.2.2 Proof of Eq. 5
  7. B Proof of Lemma 1
  8. C More details on the circuit compilation of $\mathcal{R}$
    1. C.1 Single-qudit case
      1. C.1.1 Circuit compilation
      2. C.1.2 Verification
    2. C.2 Many-qubit case
      1. C.2.1 Circuit compilation
      2. C.2.2 Verification
  9. D 2 and the proof of Theorem 2
    1. D.1 2 and its proof
    2. D.2 Proof of Theorem 2
      1. D.2.1 Unbiasedness of the estimator $\widehat{\tr_{\bar{A}}(\rho^{t})}$
      2. D.2.2 Proof of Eq. 18
    3. D.3 The estimator of moment $P_{t}$ and its variance
  10. E New estimator $\widehat{\rho^{t}}_{*}$ and its performance
  11. F Proof of Corollary 1
  12. G Details of numerical simulations

##  Appendix A Proofs of 1 and Theorem 1

###  A.1 Proof of 1

Recall that $\Pr(\mathbf{b}|V)=\bra{\mathbf{b}}V(\rho^{t})\ket{\mathbf{b}}$
and
$Q_{\mathbf{b}}=t^{-1}\sum_{i=1}^{t}\outerproduct{\mathbf{b}}{\mathbf{b}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)}$.
Hence, the first line of 1 can be proved by showing that

$\displaystyle\bra{\mathbf{b}}V(\rho^{t})\ket{\mathbf{b}}=\tr\left[S_{t}\left(\outerproduct{\mathbf{b}}{\mathbf{b}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)}\right)V(\rho)^{\otimes
t}\right]\quad\forall i=1,2,\dots,t.$ (23)

We graphically show the derivation of this relation with the diagram in Fig.
6.

The second line of 1 can be derived as follows,

$\displaystyle\tr\left(S_{t}Q_{\mathbf{b}}V(\rho)^{\otimes t}\right)$
$\displaystyle=\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}S_{t}Q_{\mathbf{b}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}$
$\displaystyle=\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}},$ (24)

where the first equality holds because
$\\{\ket{\psi_{\mathbf{x}}}\\}_{\mathbf{x}}$ forms an orthonormal basis on
$\mathcal{H}_{d}^{\otimes t}$, and the second equality holds because
$Q_{\mathbf{b}}$, $S_{t}$, and
$\outerproduct{\psi_{\mathbf{x}}}{\psi_{\mathbf{x}}}$ commute with each other
(recall that $S_{t}$ and $Q_{\mathbf{b}}$ can be diagonalized simultaneously
in the basis $\\{\ket{\psi_{\mathbf{x}}}\\}$).

Figure 6: A proof sketch of Eq. 23 for the AFRS protocol. The gray dashed
lines denote the periodic boundary condition, which is the trace operation.

###  A.2 Proof of Theorem 1

The proof of Theorem 1 is separated into two parts. Using 1, the first part in
Section A.2.1 shows that $\widehat{\rho^{t}}$ defined in Eq. 4 reproduce the
$t$-th power of the underlying state $\rho$ exactly in expectation (over the
randomness in the choice of unitary $V$, measurement outcome $\mathbf{x}$, and
the mapping strategy $\mathcal{G}_{\mathbf{x}\to\mathbf{b}}$), i.e.,
$\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\widehat{\rho^{t}}=\rho^{t}$. The
second part in Section A.2.2 proves the upper bound of
$\mathrm{Var}(\widehat{o_{t}})$ given in Eq. 5.

####  A.2.1 Unbiasedness of the estimator $\widehat{\rho^{t}}$

The expectation of $\widehat{\rho^{t}}$ is given by

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\widehat{\rho^{t}}$
$\displaystyle=\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\,{\mathrm{Re}}(f(\mathbf{x}))\cdot\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x},\mathbf{b}}\Pr(\mathbf{x},\mathbf{b}|V)\,{\mathrm{Re}}(f(\mathbf{x}))\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)\left[\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}|\mathbf{x})\,{\mathrm{Re}}(f(\mathbf{x}))\right].$
(25)

Here, the third equality follows because
$\Pr(\mathbf{x},\mathbf{b}|V)=\Pr(\mathbf{x}|V)\Pr(\mathbf{b}|\mathbf{x})$,
since the mapping result $\mathbf{b}$ only depends on the real measurement
outcome $\mathbf{x}$ according to the mapping strategy
$\mathcal{G}_{\mathbf{x}\to\mathbf{b}}$ in the main text.

Since $\Pr(\mathbf{x}|V)$ and $\Pr(\mathbf{b}|\mathbf{x})$ are real numbers,
the last term in the parentheses in Section A.2.1 can be further written as

$\displaystyle{\mathrm{Re}}\left[\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}|\mathbf{x})f(\mathbf{x})\right]={\mathrm{Re}}\left[\Pr(\mathbf{b}|V)\right]=\bra{\bf{b}}V(\rho^{t})\ket{\bf{b}},$
(26)

where the first equality follows from 1 in the main text. By plugging this
relation into Section A.2.1, we get

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\widehat{\rho^{t}}$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\bra{\bf{b}}V(\rho^{t})\ket{\bf{b}}\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)$
$\displaystyle=\mathcal{M}^{-1}\left[\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\bra{\bf{b}}V(\rho^{t})\ket{\bf{b}}V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right]=\mathcal{M}^{-1}\left[\mathcal{M}(\rho^{t})\right]=\rho^{t},$
(27)

where we used the definition of the channel $\mathcal{M}$, that is, for any
operator $X$ on $\mathcal{H}_{d}$,

$\displaystyle\mathcal{M}(X):=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\bra{\bf{b}}V(X)\ket{\bf{b}}V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V.$
(28)

In conclusion, $\widehat{\rho^{t}}$ is an unbiased estimator of $\rho^{t}$.

####  A.2.2 Proof of Eq. 5

Since
$\widehat{\rho^{t}}={\mathrm{Re}}(f(\mathbf{x}))\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)$
and $O$ are Hermitian operators, the estimator
$\widehat{o_{t}}=\tr\left(O\widehat{\rho^{t}}\right)$ is a real number, and
its variance can be written as

$\displaystyle\mathrm{Var}(\widehat{o_{t}})=\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\left[\widehat{o_{t}}^{2}\right]-\left(\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\widehat{o_{t}}\right)^{2}=\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\left[\widehat{o_{t}}^{2}\right]-\left[\tr(O\rho^{t})\right]^{2}.$
(29)

Note that the inverse of $\mathcal{M}$ in Eq. 28 is self-adjoint, i.e.,
$\tr\\!\left[X\mathcal{M}^{-1}(Y)\right]=\tr\\!\left[\mathcal{M}^{-1}(X)Y\right]$
for any pair of operators $X,Y$ on $\mathcal{H}_{d}$. So we have

$\displaystyle\widehat{o_{t}}=\tr(O\widehat{\rho^{t}})={\mathrm{Re}}(\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}})\tr\left[O\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)\right]={\mathrm{Re}}(\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}})\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}.$
(30)

It follows that

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\left[\widehat{o_{t}}^{2}\right]$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x},\mathbf{b}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}|\mathbf{x})\,\widehat{o_{t}}^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x},\mathbf{b}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}\left[{\mathrm{Re}}(\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}})\right]^{2}\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}^{2}\underbrace{\left(\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}\right)}_{(*)},$
(31)

where $(a)$ follows because
$\Pr(\mathbf{x}|V)=\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}$ and
$\Pr(\mathbf{b}|\mathbf{x})=\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}$,
and $(b)$ follows because
$\left[\mathrm{Re}(f(\mathbf{x}))\right]^{2}=\left[\mathrm{Re}(\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}})\right]^{2}\leq
1$. The term $(*)$ in Section A.2.2 can be further written as

$\displaystyle(*)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}\stackrel{{\scriptstyle(b)}}{{=}}\tr\\!\left[V(\rho)^{\otimes
t}Q_{\mathbf{b}}\right]=\bra{\mathbf{b}}V(\rho)\ket{\mathbf{b}},$ (32)

where $(a)$ holds because $Q_{\mathbf{b}}$ commutes with
$\outerproduct{\psi_{\mathbf{x}}}{\psi_{\mathbf{x}}}$, and $(b)$ holds because
$\\{\ket{\psi_{\mathbf{x}}}\\}_{\mathbf{x}}$ forms an orthonormal basis on
$\mathcal{H}_{d}^{\otimes t}$.

Equations (A.2.2) and (32) together imply that

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}\\}}\left[\widehat{o_{t}}^{2}\right]\leq\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}^{2}\bra{\mathbf{b}}V(\rho)\ket{\mathbf{b}}=\mathbb{E}\left[\operatorname{tr}(O\hat{\rho})^{2}\right],$
(33)

where $\hat{\rho}$ is the original single-copy shadow snapshot [19]. Hence,
the expectation value of $\widehat{o_{t}}^{2}$ is not larger than that of
$\left[\tr(O\hat{\rho})\right]^{2}$ in the original shadow protocol on a
single-copy of $\rho$, no matter what value $t$ takes.

By Eqs. (29) and (33), we have

$\displaystyle\mathrm{Var}(\widehat{o_{t}})$
$\displaystyle\leq\mathbb{E}\left[\operatorname{tr}(O\hat{\rho})^{2}\right]-\left[\tr(O\rho^{t})\right]^{2}$
$\displaystyle=\mathbb{V}^{\mathcal{E}}(O,\rho)+\left[\tr(O\rho)\right]^{2}-\left[\tr(O\rho^{t})\right]^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{V}^{\mathcal{E}}(O_{0},\rho)+\left[\tr(O\rho)\right]^{2}-\left[\tr(O\rho^{t})\right]^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\left\|O_{0}\right\|_{\mathrm{sh},\mathcal{E}}^{2}+\|O\|_{\infty}^{2},$
(34)

which confirms Eq. 5 in the main text. Here,
$\mathbb{V}^{\mathcal{E}}(O,\rho)$ is the variance of $\tr(O\hat{\rho})$
associated with the original shadow protocol on a single-copy of $\rho$ [19],
$(a)$ holds because shifting the operator $O$ to its traceless part $O_{0}$
does not change the variance of $\tr(O\hat{\rho})$, $(b)$ holds because
$\mathbb{V}^{\mathcal{E}}(O_{0},\rho)\leq\left\|O_{0}\right\|_{\mathrm{sh},\mathcal{E}}^{2}$
and $\left|\tr(O\rho)\right|\leq\|O\|_{\infty}$.

##  Appendix B Proof of Lemma 1

Suppose $[\mathbf{z}]$, $[\mathbf{y}]$ are two classes of $\mathbb{B}(d,t)$,
integers $k_{1}\in\\{0,1,\dots,|[\mathbf{z}]|-1\\}$ and
$k_{2}\in\\{0,1,\dots,|[\mathbf{y}]|-1\\}$. When
$[\mathbf{z}]\neq[\mathbf{y}]$, the two states
$\ket{\Psi_{[\mathbf{z}]}^{(k_{1})}}$ and
$\ket{\Psi_{[\mathbf{y}]}^{(k_{2})}}$ are orthogonal since
$\langle\tau^{r_{2}}(\mathbf{z})|\tau^{r_{1}}(\mathbf{y})\rangle=0$ for any
integers $r_{1}$ and $r_{2}$. When $[\mathbf{z}]=[\mathbf{y}]$ but $k_{1}\neq
k_{2}$, the two states $\ket{\Psi_{[\mathbf{z}]}^{(k_{1})}}$ and
$\ket{\Psi_{[\mathbf{y}]}^{(k_{2})}}$ are still orthogonal since

$\displaystyle\left\langle\Psi_{[\mathbf{z}]}^{(k_{1})}\Big{|}\Psi_{[\mathbf{y}]}^{(k_{2})}\right\rangle$
$\displaystyle=\frac{1}{|[\mathbf{z}]|}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{r(k_{2}-k_{1})\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)=0,$ (35)

where the last inequality follows because $k_{2}-k_{1}\neq
0\mod|[\mathbf{z}]|$. In addition, notice that the size of the set
$\left\\{\ket{\Psi_{[\mathbf{z}]}^{(k)}}\right\\}_{[\mathbf{z}],k}$ equals
$|\mathbb{B}(d,t)|$, and thus equals the dimension of
$\mathcal{H}_{d}^{\otimes t}$. Therefore,
$\left\\{\ket{\Psi_{[\mathbf{z}]}^{(k)}}\right\\}_{[\mathbf{z}],k}$ forms an
orthonormal basis on $\mathcal{H}_{d}^{\otimes t}$.

Recall that $S_{t}\ket{\mathbf{x}}=\ket{\tau(\mathbf{x})}$ for
$\mathbf{x}\in\mathbb{B}(d,t)$. So we have

$\displaystyle S_{t}\ket{\Psi_{[\mathbf{z}]}^{(k)}}$
$\displaystyle=\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{rk\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)S_{t}\ket{\tau^{r}(\mathbf{z})}$
$\displaystyle=\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{rk\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\ket{\tau^{r+1}(\mathbf{z})}$
$\displaystyle=\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r^{\prime}=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{(r^{\prime}-1)k\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\ket{\tau^{r^{\prime}}(\mathbf{z})}$
$\displaystyle=\exp\left(-\frac{k\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r^{\prime}=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{r^{\prime}k\cdot
2\pi i}{|[\mathbf{z}]|}\right)\ket{\tau^{r^{\prime}}(\mathbf{z})}$
$\displaystyle=\exp\left(-\frac{k\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\ket{\Psi_{[\mathbf{z}]}^{(k)}}.$ (36)

Hence, $\ket{\Psi_{[\mathbf{z}]}^{(k)}}$ is an eigenstate of $S_{t}$ with
eigenvalue $\exp\left(-\frac{k\cdot 2\pi\rm{i}}{|[\mathbf{z}]|}\right)$.

For $\mathbf{x}\in\mathbb{B}(d,t)$, we have

$\displaystyle\left(\outerproduct{\mathbf{b}}{\mathbf{b}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)}\right)\ket{\mathbf{x}}=\begin{cases}\ket{\mathbf{x}}&\text{the
$i$-th number of $\mathbf{x}$ equals $\mathbf{b}$},\\\
0&\text{otherwise}.\end{cases}$ (37)

Consequently, we have

$\displaystyle
Q_{\mathbf{b}}\ket{\mathbf{x}}=\frac{1}{t}\sum_{i}\left(\outerproduct{\mathbf{b}}{\mathbf{b}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)}\right)\ket{\mathbf{x}}=\frac{\text{the
number of $\mathbf{b}$ in the string $\mathbf{x}$}}{t}\ket{\mathbf{x}},$ (38)

and

$\displaystyle Q_{\mathbf{b}}\ket{\Psi_{[\mathbf{z}]}^{(k)}}$
$\displaystyle=\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{rk\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)Q_{\mathbf{b}}\ket{\tau^{r}(\mathbf{z})}$
$\displaystyle=\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{rk\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\frac{\text{the number of $\mathbf{b}$ in
the string $\tau^{r}(\mathbf{z})$}}{t}\ket{\tau^{r}(\mathbf{z})}$
$\displaystyle=\frac{\text{the number of $\mathbf{b}$ in the string
$\mathbf{z}$}}{t}\times\frac{1}{\sqrt{|[\mathbf{z}]|}}\sum_{r=0}^{|[\mathbf{z}]|-1}\exp\left(\frac{rk\cdot
2\pi\rm{i}}{|[\mathbf{z}]|}\right)\ket{\tau^{r}(\mathbf{z})}$
$\displaystyle=\frac{\text{the number of $\mathbf{b}$ in the string
$\mathbf{z}$}}{t}\ket{\Psi_{[\mathbf{z}]}^{(k)}},$ (39)

where the third equality holds because the strings $\mathbf{z}$ and
$\tau^{r}(\mathbf{z})$ contains the same number of $\mathbf{b}$ for any
integer $r$. Hence, $\ket{\Psi_{[\mathbf{z}]}^{(k)}}$ is an eigenstate of
$Q_{\mathbf{b}}$ with eigenvalue $t^{-1}\times($the number of $\mathbf{b}$ in
the string $\mathbf{z})$. This completes the proof.

##  Appendix C More details on the circuit compilation of $\mathcal{R}$

###  C.1 Single-qudit case

####  C.1.1 Circuit compilation

One can follow the route of the single-qubit case in the main text to find the
circuit compilation of $\mathcal{R}$ in the general single-qudit cases, where
$\mathcal{R}$ acting on $\mathcal{H}_{d}^{\otimes 2}$ has a matrix form given
in Eq. 8. Similar as the single-qubit case, the key idea is to use measurement
on the second replica to effectively project (partially decoherence) the input
into one of the subspaces $\mathcal{V}_{c}={\rm
span}\left(\\{\ket{ab}\\}_{a+b=c}\right)$ for $c=0,1,\dots,d-2$. Note that
each $\mathcal{V}_{c}$ can be written as a direct sum of several 1- or
2-dimensional subspace $\mathcal{V}_{a,b}={\rm
span}\left(\ket{ab},\ket{ba}\right)$, i.e.,

$\displaystyle\mathcal{V}_{c}=\bigoplus_{\begin{subarray}{c}a,b\in\\{0,1,\dots,d-1\\}\\\
a+b=c\end{subarray}}\mathcal{V}_{a,b}.$ (40)

Hence, after the measurement, one can then create the desired superposition on
each $\mathcal{V}_{a,b}$ by performing some unitary transformation on the
first replica.

Figure 7: Circuit compilation of $\mathcal{R}$ in the single-qudit case.

The procedure is illustrated in Fig. 7. Different from the single-qubit case,
here we expand the space dimension of the second replica from $d$ to $2d$, and
perform a $CX$ operation from $a$ to $b$, which transforms the second replica
from $b$ to $a+b$. Next, we measure the second replica in the computational
basis $\\{|i\rangle\\}_{i=0}^{2d-1}$, and get the measurement result $c=a+b$.
Define $q_{1}:=\max\\{0,c-d+1\\}$ and $q_{2}:=\min\\{c,d-1\\}$. The form of
the unitary $U_{c}$ controlled by the measurement result $c$ depends on the
parity of $c$. For odd $c$, we have

$\displaystyle U_{c}=$
$\displaystyle\sum_{j=q_{1}}^{\frac{c-1}{2}}\left(\frac{|j\rangle+|c-j\rangle}{\sqrt{2}}\right)\langle
j|+\sum_{j=\frac{c+1}{2}}^{q_{2}}\left(\frac{|c-j\rangle-|j\rangle}{\sqrt{2}}\right)\langle
j|+\sum_{j<q_{1},j>q_{2}}\outerproduct{j}{j}.$ (41)

For even $c$, we have

$\displaystyle U_{c}=$
$\displaystyle\sum_{j=q_{1}}^{\frac{c}{2}-1}\left(\frac{|j\rangle+|c-j\rangle}{\sqrt{2}}\right)\langle
j|+\outerproduct{\frac{c}{2}}{\frac{c}{2}}+\sum_{j=\frac{c}{2}+1}^{q_{2}}\left(\frac{|c-j\rangle-|j\rangle}{\sqrt{2}}\right)\langle
j|+\sum_{j<q_{1},j>q_{2}}\outerproduct{j}{j}.$ (42)

At last, we measure the first replica in the computational basis
$\\{|i\rangle\\}_{i=0}^{d-1}$, and denote the measurement result as $x_{1}$.
By classical post-processing, we let $x_{1}$ and $x_{2}=c-x_{1}$ be the final
output of the first and second replica, respectively. Note that $U_{c}$ only
introduces only some superposition on selected pairs of computational bases of
a single-qudit, and we remark such operations could be realized, for example,
in optical systems by path encoding.

####  C.1.2 Verification

To verify that our circuit indeed implements the desired unitary $\mathcal{R}$
and the subsequent measurement, let us consider the evolution of the 2-qudit
joint system. The initial state of the joint system can be expressed as
$\ket{\Phi_{0}}=\sum_{a,b=0}^{d-1}t_{a,b}\ket{ab}$. After the $CX$ operation,
the joint state evolves into
$\ket{\Phi_{1}}=\sum_{a,b=0}^{d-1}t_{a,b}\ket{a}\ket{a+b}$. After measuring
the second qudit, depending on the measurement outcome $c$, the reduced state
on the first replica reads

$\displaystyle\ket{\Phi_{2}}=\frac{(\mathbb{I}\otimes\bra{c})\ket{\Phi_{1}}}{\sqrt{\Pr_{2}(c)}}=\frac{1}{\sqrt{\Pr_{2}(c)}}\sum_{j=q_{1}}^{q_{2}}t_{j,c-j}\ket{j},$
(43)

where

$\displaystyle{\Pr}_{2}(c)=\bra{\Phi_{1}}\left(\mathbb{I}\otimes\outerproduct{c}{c}\right)\ket{\Phi_{1}}=\sum_{j=q_{1}}^{q_{2}}\left|t_{j,c-j}\right|^{2}$
(44)

is the probability of obtaining outcome $c$. Next, we consider two cases
depending on the parity of $c$.

  * 1.

The measurement outcome $c$ is odd.
After the unitary evolution $U_{c}$, the state on the first replica reads

$\displaystyle\ket{\Phi_{3}}=U_{c}\ket{\Phi_{2}}=\frac{1}{\sqrt{\Pr_{2}(c)}}\left(\sum_{j=q_{1}}^{\frac{c-1}{2}}t_{j,c-j}\frac{|j\rangle+|c-j\rangle}{\sqrt{2}}+\sum_{j=\frac{c+1}{2}}^{q_{2}}t_{j,c-j}\frac{|c-j\rangle-|j\rangle}{\sqrt{2}}\right).$
(45)

If we measure the first replica at this point, then the conditional
probability of obtaining outcome $x_{1}$ is given by

$\displaystyle{\Pr}_{1}(x_{1}|c)=\left|\langle
x_{1}\\!\ket{\Phi_{3}}\right|^{2}=\frac{1}{\Pr_{2}(c)}\times\left\\{\begin{array}[]{lll}0&\
\text{when $x_{1}<q_{1}$ or $x_{1}>q_{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{c-x_{1},x_{1}}+t_{x_{1},c-x_{1}}\right|^{2}&\ \text{when
$q_{1}\leq x_{1}\leq\frac{c-1}{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{c-x_{1},x_{1}}-t_{x_{1},c-x_{1}}\right|^{2}&\ \text{when
$\frac{c+1}{2}\leq x_{1}\leq q_{2}$}.\end{array}\right.$ (49)

  * 2.

The measurement outcome $c$ is even.
After the unitary evolution $U_{c}$, the state on the first replica reads

$\displaystyle\ket{\Phi_{3}}=U_{c}\ket{\Phi_{2}}=\frac{1}{\sqrt{\Pr_{2}(c)}}\left(\sum_{j=q_{1}}^{\frac{c}{2}-1}t_{j,c-j}\frac{|j\rangle+|c-j\rangle}{\sqrt{2}}+t_{\frac{c}{2},\frac{c}{2}}\ket{\frac{c}{2}}+\sum_{j=\frac{c}{2}+1}^{q_{2}}t_{j,c-j}\frac{|c-j\rangle-|j\rangle}{\sqrt{2}}\right).$
(50)

If we measure the first replica at this point, then the conditional
probability of obtaining outcome $x_{1}$ is given by

$\displaystyle{\Pr}_{1}(x_{1}|c)=\left|\langle
x_{1}\\!\ket{\Phi_{3}}\right|^{2}=\frac{1}{\Pr_{2}(c)}\times\left\\{\begin{array}[]{lll}0&\
\text{when $x_{1}<q_{1}$ or $x_{1}>q_{2}$},\\\\[2.15277pt]
|t_{\frac{c}{2},\frac{c}{2}}|^{2}&\ \text{when
$x_{1}=\frac{c}{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{c-x_{1},x_{1}}+t_{x_{1},c-x_{1}}\right|^{2}&\ \text{when
$q_{1}\leq x_{1}\leq\frac{c}{2}-1$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{c-x_{1},x_{1}}-t_{x_{1},c-x_{1}}\right|^{2}&\ \text{when
$\frac{c}{2}+1\leq x_{1}\leq q_{2}$}.\end{array}\right.$ (55)

Recall that the final result for the second qudit is $x_{2}=c-x_{1}$.
According to Eqs. (49) and (55), for our circuit, the overall probability of
obtaining outcomes $x_{1}$ on the first qudit and $x_{2}$ on the second qudit
reads

$\displaystyle\Pr(x_{1},x_{2})$
$\displaystyle={\Pr}_{1}(x_{1}|c=x_{1}+x_{2})\,{\Pr}_{2}(c=x_{1}+x_{2})=\left\\{\begin{array}[]{lll}|t_{x_{1},x_{1}}|^{2}&\
\text{when $x_{1}=x_{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{x_{2},x_{1}}+t_{x_{1},x_{2}}\right|^{2}&\ \text{when
$x_{1}<x_{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{x_{2},x_{1}}-t_{x_{1},x_{2}}\right|^{2}&\ \text{when
$x_{1}>x_{2}$}.\end{array}\right.$ (59)

On the other hand, if the initial state $\ket{\Phi_{0}}$ is evolved under
$\mathcal{R}$, then it becomes

$\displaystyle\ket{\Phi_{1}^{\prime}}=\mathcal{R}\ket{\Phi_{0}}=\sum_{a=0}^{d-1}t_{a,a}\ket{aa}+\sum_{a<b}t_{a,b}\frac{\ket{ba}+\ket{ab}}{\sqrt{2}}+\sum_{a>b}t_{a,b}\frac{\ket{ba}-\ket{ab}}{\sqrt{2}}.$
(60)

If we measure this state in the computational basis, then the probability of
obtaining outcomes $x_{1}$ on the first qudit and $x_{2}$ on the second qudit
reads

$\displaystyle\left|\langle
x_{1},x_{2}\\!\ket{\Phi_{1}^{\prime}}\right|^{2}=\left\\{\begin{array}[]{lll}|t_{x_{1},x_{1}}|^{2}&\
\text{when $x_{1}=x_{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{x_{2},x_{1}}+t_{x_{1},x_{2}}\right|^{2}&\ \text{when
$x_{1}<x_{2}$},\\\\[2.15277pt]
\frac{1}{2}\left|t_{x_{2},x_{1}}-t_{x_{1},x_{2}}\right|^{2}&\ \text{when
$x_{1}>x_{2}$}.\end{array}\right.$ (64)

This probability exactly equals the outcome probability of our circuit shown
in Eq. 59. Due to linearity, this result also holds when the initial state of
the 2-qudit system is mixed. Therefore, our circuit in Fig. 7 can indeed
realize the desired unitary $\mathcal{R}$ and the subsequent measurement.

###  C.2 Many-qubit case

####  C.2.1 Circuit compilation

Here, we give more explanations on our circuit compilation of $\mathcal{R}$ in
the many-qubit case. We use $\mathbf{a}$ and $\mathbf{b}$ to label the
$n$-qubit basis for the first and second replicas, respectively. As mentioned
in the main text, the action of $\mathcal{R}$ on
$|\mathbf{a}\rangle|\mathbf{b}\rangle$ is as follows: if $a^{i}=b^{i}$,
$\mathcal{R}$ does nothing; and for the rest of qubits in $[l]$ with
$a^{i}\neq b^{i}$, $R$ creates a GHZ-like ‘large’ superposition in the
subspace $\mathcal{V}_{\mathbf{a}^{[l]}}={\rm
span}\\!\left(\ket{\mathbf{a}^{[l]}}\ket{\overline{\mathbf{a}^{[l]}}},\ket{\overline{\mathbf{a}^{[l]}}}\ket{\mathbf{a}^{[l]}}\right)$.

Following the circuit compilation in the single-qubit case, our circuit for
$\mathcal{R}$ in the many-qubit case is designed as follows; see Fig. 8 (a)
for an illustration. Still, we use qubits of the second replica to check the
parity information between the two replicas. This is again achieved by first
applying a $CX_{a^{i}\rightarrow b^{i}}$ for each qubit pair $(a^{i},b^{i})$,
and then measuring the second qubit in each pair in the computational basis.
If the measurement result is $c^{i}=0$, then we know that $a^{i}=b^{i}$.
Otherwise, if the measurement result is $c^{i}=1$, we conclude that $a^{i}\neq
b^{i}$. Let $[l]:=\\{i\,|\,c^{i}=1\\}$. Notice that the measurement
effectively projects (partially decoherence) the qubit pairs in $[l]$ into the
following subspace

$\displaystyle\mathcal{V}_{[l]}=\bigoplus_{\mathbf{a}^{[l]}\in\\{0,1\\}^{l}}\mathcal{V}_{\mathbf{a}^{[l]}}.$
(65)

Hence, after the measurement, one can then create the desired superposition on
each $\mathcal{V}_{\mathbf{a}^{[l]}}$ by performing a Hadamard-like
transformation $U_{[l]}^{\prime}$ on all qubits $\\{a^{i}\\}_{i\in[l]}$. The
action of $U_{[l]}^{\prime}$ is given as

$\displaystyle
U_{[l]}^{\prime}\ket{\mathbf{a}^{[l]}}=\begin{cases}\frac{1}{\sqrt{2}}\left(\ket{\mathbf{a}^{[l]}}+\ket{\overline{\mathbf{a}^{[l]}}}\right)&\quad\mathbf{a}^{[l]}<\overline{\mathbf{a}^{[l]}},\\\\[4.30554pt]
\frac{1}{\sqrt{2}}\left(\ket{\overline{\mathbf{a}^{[l]}}}-\ket{\mathbf{a}^{[l]}}\right)&\quad\mathbf{a}^{[l]}>\overline{\mathbf{a}^{[l]}}.\end{cases}$
(66)

At last, we measure all $n$ qubits of the first replica in the computational
basis. As shown in Fig. 8 (a), for the $i$-th qubit, we denote the measurement
result as $x_{1}^{i}$. In classical post-processing, we let $x_{1}^{i}$ and
$x_{2}^{i}=c^{i}\oplus x_{1}^{i}$ be the final output of the $i$-th qubit
pair.

Figure 8: (a) Circuit compilation of $\mathcal{R}$ in the many-qubit case.
Here, we take $n=4$, and assume that the parity check $c^{1}=0$, and $c^{i}=1$
for $i\in\\{2,3,4\\}$. (b) Gate synthesis of $U_{[l]}^{\prime}$.

The remaining task is to decompose the unitary $U_{[l]}^{\prime}$ into
elementary quantum gates. To this end, we consider the subspace ${\rm
span}\\!\left(\ket{\mathbf{a}^{[l]}},\ket{\overline{\mathbf{a}^{[l]}}}\right)$
as a code space, where $\ket{0_{L}}=\ket{\mathbf{a}^{[l]}}$ and
$\ket{1_{L}}=\ket{\overline{\mathbf{a}^{[l]}}}$ are the logical 0 and 1
states, respectively (here we assume
$\mathbf{a}^{[l]}<\overline{\mathbf{a}^{[l]}}$ without lose of generality).
Recall that the condition $\mathbf{a}^{[l]}<\overline{\mathbf{a}^{[l]}}$ means
that the first bit in the string $\mathbf{a}^{[l]}$
$\left(\overline{\mathbf{a}^{[l]}}\right)$ is 0 (1). Using this fact, one can
easily verify that the operations $X^{\otimes l}$ and $Z_{1}$ act as logical
$X$ and $Z$ gates on the code space, where $Z_{1}$ denotes the Pauli-$Z$
operator on the first qubit. According to Eq. 66, the unitary
$U_{[l]}^{\prime}$ acts as a logical Hadamard operation on the code space.
Therefore, it can be chosen as the following form,

$\displaystyle U_{[l]}^{\prime}=\frac{1}{\sqrt{2}}\left(X^{\otimes
l}+Z_{1}\right),$ (67)

which is a Clifford operator. Now it suffices to find another Clifford unitary
$\tilde{U}_{[l]}$ such that
$\tilde{U}_{[l]}H_{1}\tilde{U}_{[l]}^{{\dagger}}=U_{[l]}^{\prime}$, where
$H_{1}=(X_{1}+Z_{1})/\sqrt{2}$ is the Hadamard gate on the first qubit. One
can check that the goal can be achieved by constructing $\tilde{U}_{[l]}$ from
a series of $CX$ gates, that is,

$\displaystyle\begin{aligned} \tilde{U}_{[l]}=CX_{l-1\rightarrow
l}\,CX_{l-2\rightarrow l-1}\cdots CX_{1\rightarrow 2}.\end{aligned}$ (68)

With this construction, the gate synthesis of $U_{[l]}^{\prime}$ is shown in
Fig. 8 (b).

By plugging the gate synthesis of $U_{[l]}^{\prime}$ into the circuit of Fig.
8 (a), we are very close to our final circuit. It worth noting that the $CX$
gates after $H_{1}$ and before the final measurement can be substituted by
classical post-processing. The final circuit after this simplification is
shown in Fig. 3 (c) of the main text.

####  C.2.2 Verification

To verify that our circuit constructed in the previous section indeed
implements the desired unitary $\mathcal{R}$ and the subsequent measurement,
let us consider the evolution of the $2n$-qubit joint system. The initial
state of the joint system can be expressed as
$\ket{\Phi_{0}}=\sum_{\textbf{a},\textbf{b}\in\\{0,1\\}^{n}}t_{\textbf{a},\textbf{b}}\ket{\textbf{a}}_{1}\ket{\textbf{b}}_{2}$.
After the $CX_{a_{i}\rightarrow b_{i}}$ operations, the joint state evolves
into
$\ket{\Phi_{1}}=\sum_{\textbf{a},\textbf{b}\in\\{0,1\\}^{n}}t_{\textbf{a},\textbf{b}}\ket{\textbf{a}}_{1}\ket{\textbf{a}\oplus\textbf{b}}_{2}$.
After measuring the $n$ qubits of the second replica, depending on the
measurement outcome $\mathbf{c}\in\\{0,1\\}^{n}$, the reduced state on the
first replica reads

$\displaystyle\ket{\Phi_{2}^{\mathbf{c}}}=\frac{(\mathbb{I}\otimes\bra{\mathbf{c}})\ket{\Phi_{1}}}{\sqrt{\Pr_{2}(\mathbf{c})}}=\frac{1}{\sqrt{\Pr_{2}(\mathbf{c})}}\sum_{\textbf{a},\textbf{b}\in\\{0,1\\}^{n}}t_{\textbf{a},\textbf{b}}\,\delta_{\textbf{c},\textbf{a}\oplus\textbf{b}}\ket{\textbf{a}}=\frac{1}{\sqrt{\Pr_{2}(\mathbf{c})}}\sum_{\textbf{a}\in\\{0,1\\}^{n}}t_{\textbf{a},\textbf{a}\oplus\textbf{c}}\ket{\textbf{a}},$
(69)

where

$\displaystyle{\Pr}_{2}(\mathbf{c})=\bra{\Phi_{1}}\left(\mathbb{I}\otimes\outerproduct{\mathbf{c}}{\mathbf{c}}\right)\ket{\Phi_{1}}=\sum_{\textbf{a}\in\\{0,1\\}^{n}}\left|t_{\textbf{a},\textbf{a}\oplus\textbf{c}}\right|^{2}$
(70)

is the probability of obtaining outcome $\mathbf{c}$. Next, we consider two
cases depending on the outcome $\mathbf{c}$.

  * 1.

The outcome $\mathbf{c}=\mathbf{0}$.
In this case, we then measure $\ket{\Phi_{2}^{\mathbf{c}}}$ directly in the
computational basis. The conditional probability of obtaining outcome
$\mathbf{x}\in\\{0,1\\}^{n}$ is given by

$\displaystyle{\Pr}_{1}(\mathbf{x}|\mathbf{c}=\mathbf{0})=\left|\langle\mathbf{x}\\!\ket{\Phi_{2}^{\mathbf{c}}}\right|^{2}=\left|\frac{1}{\sqrt{\Pr_{2}(\mathbf{c})}}\sum_{\textbf{a}\in\\{0,1\\}^{n}}t_{\textbf{a},\textbf{a}}\delta_{\mathbf{a},\mathbf{x}}\right|^{2}=\frac{\left|t_{\mathbf{x},\mathbf{x}}\right|^{2}}{\Pr_{2}(\mathbf{c})}.$
(71)

  * 2.

The outcome $\mathbf{c}\neq\mathbf{0}$.
Let $[k]:=\\{i\,|\,c^{i}=0\\}$ and $[l]:=\\{i\,|\,c^{i}=1\\}$ with $k+l=n$.
For any bit string $\textbf{a}\in\\{0,1\\}^{n}$, we write
$\mathbf{a}=\mathbf{a}^{[l]}\mathbf{a}^{[k]}$. In this case, the unitary
$U_{[l]}^{\prime}$ is applied on the qubits of $\ket{\Phi_{2}^{\mathbf{c}}}$
in $[l]$, which evolves the state into

$\displaystyle\ket{\Phi_{3}^{\mathbf{c}}}$
$\displaystyle=\left(\mathbb{I}_{[k]}\otimes
U_{[l]}^{\prime}\right)\ket{\Phi_{2}^{\mathbf{c}}}$
$\displaystyle=\frac{1}{\sqrt{\Pr_{2}(\mathbf{c})}}\sum_{\textbf{a}\in\\{0,1\\}^{n}}t_{\textbf{a},\textbf{a}\oplus\textbf{c}}\ket{\textbf{a}^{[k]}}\otimes\left(U_{[l]}^{\prime}\ket{\textbf{a}^{[l]}}\right)$
$\displaystyle=\frac{1}{\sqrt{\Pr_{2}(\mathbf{c})}}\sum_{\mathbf{a}^{[k]}\in\\{0,1\\}^{n-l}}\ket{\textbf{a}^{[k]}}\otimes\left[\sum_{\mathbf{a}^{[l]}|\,\mathbf{a}^{[l]}<\overline{\mathbf{a}^{[l]}}}t_{\textbf{a},\textbf{a}\oplus\textbf{c}}\,\frac{\ket{\mathbf{a}^{[l]}}+\ket{\overline{\mathbf{a}^{[l]}}}}{\sqrt{2}}+\sum_{\mathbf{a}^{[l]}|\,\mathbf{a}^{[l]}>\overline{\mathbf{a}^{[l]}}}t_{\textbf{a},\textbf{a}\oplus\textbf{c}}\,\frac{\ket{\overline{\mathbf{a}^{[l]}}}-\ket{\mathbf{a}^{[l]}}}{\sqrt{2}}\right].$
(72)

If we measure the first replica at this point, then the conditional
probability of obtaining outcome $\mathbf{x}\in\\{0,1\\}^{n}$ is given by

$\displaystyle{\Pr}_{1}(\mathbf{x}|\mathbf{c})=\left|\langle\mathbf{x}\\!\ket{\Phi_{3}^{\mathbf{c}}}\right|^{2}$
$\displaystyle=\frac{1}{\Pr_{2}(\mathbf{c})}\times\left\\{\begin{array}[]{lll}\frac{1}{2}\left|t_{\overline{\mathbf{x}^{[l]}}\mathbf{x}^{[k]},\mathbf{x}^{[l]}\mathbf{x}^{[k]}}+t_{\mathbf{x}^{[l]}\mathbf{x}^{[k]},\overline{\mathbf{x}^{[l]}}\mathbf{x}^{[k]}}\right|^{2}&\
\text{when $\mathbf{x}^{[l]}<\overline{\mathbf{x}^{[l]}}$},\\\\[5.16663pt]
\frac{1}{2}\left|t_{\overline{\mathbf{x}^{[l]}}\mathbf{x}^{[k]},\mathbf{x}^{[l]}\mathbf{x}^{[k]}}-t_{\mathbf{x}^{[l]}\mathbf{x}^{[k]},\overline{\mathbf{x}^{[l]}}\mathbf{x}^{[k]}}\right|^{2}&\
\text{when $\mathbf{x}^{[l]}>\overline{\mathbf{x}^{[l]}}$},\end{array}\right.$
(75)
$\displaystyle=\frac{1}{\Pr_{2}(\mathbf{c})}\times\left\\{\begin{array}[]{lll}\frac{1}{2}\left|t_{\mathbf{x}\oplus\mathbf{c},\mathbf{x}}+t_{\mathbf{x},\mathbf{x}\oplus\mathbf{c}}\right|^{2}&\
\text{when $\mathbf{x}<\mathbf{x}\oplus\mathbf{c}$},\\\\[3.01385pt]
\frac{1}{2}\left|t_{\mathbf{x}\oplus\mathbf{c},\mathbf{x}}-t_{\mathbf{x},\mathbf{x}\oplus\mathbf{c}}\right|^{2}&\
\text{when $\mathbf{x}>\mathbf{x}\oplus\mathbf{c}$},\end{array}\right.$ (78)

where we have used the relation
$\mathbf{x}\oplus\mathbf{c}=\overline{\mathbf{x}^{[l]}}\mathbf{x}^{[k]}$.

Recall that the final output result for the second replica is
$\mathbf{x}\oplus\mathbf{c}$ if the measurement outcome of the first replica
is $\mathbf{x}$. According to Eqs. (71) and (75), for our circuit, the overall
probability of obtaining outcomes $\mathbf{x}_{1}\in\\{0,1\\}^{n}$ on the
first replica and $\mathbf{x}_{2}\in\\{0,1\\}^{n}$ on the second replica reads

$\displaystyle\Pr(\mathbf{x}_{1},\mathbf{x}_{2})$
$\displaystyle={\Pr}_{1}(\mathbf{x}_{1}|\mathbf{c}=\mathbf{x}_{1}\oplus\mathbf{x}_{2})\,{\Pr}_{2}(\mathbf{c}=\mathbf{x}_{1}\oplus\mathbf{x}_{2})=\left\\{\begin{array}[]{lll}\left|t_{\mathbf{x}_{1},\mathbf{x}_{1}}\right|^{2}.&\
\text{when $\mathbf{x}_{1}=\mathbf{x}_{2}$},\\\\[3.44444pt]
\frac{1}{2}\left|t_{\mathbf{x}_{2},\mathbf{x}_{1}}+t_{\mathbf{x}_{1},\mathbf{x}_{2}}\right|^{2}&\
\text{when $\mathbf{x}_{1}<\mathbf{x}_{2}$},\\\\[3.44444pt]
\frac{1}{2}\left|t_{\mathbf{x}_{2},\mathbf{x}_{1}}-t_{\mathbf{x}_{1},\mathbf{x}_{2}}\right|^{2}&\
\text{when $\mathbf{x}_{1}>\mathbf{x}_{2}$}.\end{array}\right.$ (82)

On the other hand, if the initial state $\ket{\Phi_{0}}$ is evolved under
$\mathcal{R}$, then it becomes

$\displaystyle\ket{\Phi_{1}^{\prime}}=\mathcal{R}\ket{\Phi_{0}}=\sum_{\mathbf{a}\in\\{0,1\\}^{n}}t_{\mathbf{a},\mathbf{a}}\ket{\mathbf{a}}_{1}\ket{\mathbf{a}}_{2}+\sum_{\mathbf{a},\mathbf{b}|\,\mathbf{a}<\mathbf{b}}t_{\mathbf{a},\mathbf{b}}\frac{\ket{\mathbf{b}}_{1}\ket{\mathbf{a}}_{2}+\ket{\mathbf{a}}_{1}\ket{\mathbf{b}}_{2}}{\sqrt{2}}+\sum_{\mathbf{a},\mathbf{b}|\,\mathbf{a}>\mathbf{b}}t_{\mathbf{a},\mathbf{b}}\frac{\ket{\mathbf{b}}_{1}\ket{\mathbf{a}}_{2}-\ket{\mathbf{a}}_{1}\ket{\mathbf{b}}_{2}}{\sqrt{2}}.$
(83)

If we measure this state in the computational basis, then the probability of
obtaining outcomes $\mathbf{x}_{1}\in\\{0,1\\}^{n}$ on the first replica and
$\mathbf{x}_{2}\in\\{0,1\\}^{n}$ on the second replica reads

$\displaystyle\left|(\bra{\mathbf{x}_{1}}\otimes\bra{\mathbf{x}_{2}})\\!\ket{\Phi_{1}^{\prime}}\right|^{2}=\left\\{\begin{array}[]{lll}\left|t_{\mathbf{x}_{1},\mathbf{x}_{1}}\right|^{2}.&\
\text{when $\mathbf{x}_{1}=\mathbf{x}_{2}$},\\\\[3.44444pt]
\frac{1}{2}\left|t_{\mathbf{x}_{2},\mathbf{x}_{1}}+t_{\mathbf{x}_{1},\mathbf{x}_{2}}\right|^{2}&\
\text{when $\mathbf{x}_{1}<\mathbf{x}_{2}$},\\\\[3.44444pt]
\frac{1}{2}\left|t_{\mathbf{x}_{2},\mathbf{x}_{1}}-t_{\mathbf{x}_{1},\mathbf{x}_{2}}\right|^{2}&\
\text{when $\mathbf{x}_{1}>\mathbf{x}_{2}$}.\end{array}\right.$ (87)

This probability exactly equals the outcome probability of our circuit shown
in Eq. 82. Due to linearity, this result also hold when the initial state of
the $2n$-qubit system is mixed. Therefore, our circuit can indeed realize the
desired unitary $\mathcal{R}$ and the subsequent measurement in the many qubit
case.

##  Appendix D 2 and the proof of Theorem 2

The goal of this appendix is to derive Theorem 2 in the main text. In
preparation for presenting the formal proof in Section D.2, we first provide
an observation in Section D.1, which clarifies the relation between the fake
conditional probability $\Pr(\mathbf{b}^{A}|V_{A})$, function
$f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})$, mapping probability
$\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})$, and the probability $\Pr(\mathbf{x}|V)$.
This observation is the counterpart of 1 in the main text. In the following,
we denote by
$\ket{\phi_{\mathbf{x}}}:=\ket{\psi_{\mathbf{x}^{A}}}\otimes\ket{\psi_{\mathbf{x}^{\bar{A}}}}$
for simplicity, where $\mathbf{x}=\mathbf{x}^{A}\mathbf{x}^{\bar{A}}$.

###  D.1 2 and its proof

######  Observation 2.

The fake and marginal conditional probability $\Pr(\mathbf{b}^{A}|V_{A})$ can
be written as

$\displaystyle\Pr(\mathbf{b}^{A}|V_{A})$
$\displaystyle=\tr\left[S_{t}\left(Q_{\mathbf{b}^{A}}\otimes{\mathbb{I}}_{\bar{A}}^{\otimes
t}\right)V(\rho)^{\otimes t}\right]$ (88)
$\displaystyle=\sum_{\mathbf{x}}\bra{\phi_{\mathbf{x}}}S_{t}\ket{\phi_{\mathbf{x}}}\bra{\phi_{\mathbf{x}}}\left(Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes
t}\right)\ket{\phi_{\mathbf{x}}}\bra{\phi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\phi_{\mathbf{x}}}$
$\displaystyle=\sum_{\mathbf{x}}f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})\Pr(\mathbf{x}|V).$

Here, $V=V_{A}\otimes V_{\bar{A}}$, and the function
$f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})$ and the mapping probability
$\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})$ are defined in Eqs. (15) and (16),
respectively.

###### Proof.

Recall that
$\Pr(\mathbf{b}^{A}|V_{A})=\bra{\mathbf{b}^{A}}V_{A}(\tr_{\bar{A}}(\rho^{t}))\ket{\mathbf{b}^{A}}$
and
$Q_{\mathbf{b}^{A}}=t^{-1}\sum_{i}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}_{i}\otimes\mathbb{I}^{\otimes(t-1)}_{A}$.
So the first line of Eq. 88 can be proved by showing that

$\bra{\mathbf{b}^{A}}V_{A}(\tr_{\bar{A}}(\rho^{t}))\ket{\mathbf{b}^{A}}=\tr\left(S_{t}\left[\left(\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}_{i}\otimes{\mathbb{I}}_{A}^{\otimes(t-1)}\right)\otimes{\mathbb{I}}_{\bar{A}}^{\otimes
t}\right]V(\rho)^{\otimes t}\right)\quad\forall i=1,2,\cdots,t.$ (89)

And we graphically show the derivation of this relation in Fig. 9.

Figure 9: A proof sketch of Eq. 89 for local-AFRS protocol. Here, the yellow
lines represent the subsystem $\bar{A}$ of $\rho$, and the gray dashed lines
denote the trace operation. It indicates in the right sketch that
$V_{\bar{A}}$ is not essential as
$V_{\bar{A}}V_{\bar{A}}^{\dagger}={\mathbb{I}}_{\bar{A}}$. Thus, in practise,
we can set $V_{\bar{A}}={\mathbb{I}}_{\bar{A}}$ for simplicity.

Recall that $S_{t}^{A}$ and $Q_{\mathbf{b}^{A}}$ can be diagonalized
simultaneously under the basis $\left\\{\ket{\psi_{\mathbf{x}^{A}}}\right\\}$,
and $S_{t}^{\bar{A}}$ can be diagonalized in the basis
$\left\\{\ket{\psi_{\mathbf{x}^{\bar{A}}}}\right\\}$. Hence, the operators
$\outerproduct{\phi_{\mathbf{x}}}{\phi_{\mathbf{x}}}=\outerproduct{\psi_{\mathbf{x}^{A}}}{\psi_{\mathbf{x}^{A}}}\otimes\outerproduct{\psi_{\mathbf{x}^{\bar{A}}}}{\psi_{\mathbf{x}^{\bar{A}}}}$,
$S_{t}=S_{t}^{A}\otimes S_{t}^{\bar{A}}$, and
$Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes t}$ commute with each
other. Using this fact, the second line of Eq. 88 can be derived as

$\displaystyle\tr\left[S_{t}\left(Q_{\mathbf{b}^{A}}\otimes{\mathbb{I}}_{\bar{A}}^{\otimes
t}\right)V(\rho)^{\otimes t}\right]$
$\displaystyle=\sum_{\mathbf{x}}\bra{\phi_{\mathbf{x}}}S_{t}\left(Q_{\mathbf{b}^{A}}\otimes{\mathbb{I}}_{\bar{A}}^{\otimes
t}\right)V(\rho)^{\otimes t}\ket{\phi_{\mathbf{x}}}$
$\displaystyle=\sum_{\mathbf{x}}\bra{\phi_{\mathbf{x}}}S_{t}\ket{\phi_{\mathbf{x}}}\bra{\phi_{\mathbf{x}}}\left(Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes
t}\right)\ket{\phi_{\mathbf{x}}}\bra{\phi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\phi_{\mathbf{x}}}.$ (90)

Finally, the third line of Eq. 88 follows because

$\displaystyle\bra{\phi_{\mathbf{x}}}\left(Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes
t}\right)\ket{\phi_{\mathbf{x}}}=\bra{\psi_{\mathbf{x}^{A}}}Q_{\mathbf{b}^{A}}\ket{\psi_{\mathbf{x}^{A}}}\bra{\psi_{\mathbf{x}^{\bar{A}}}}\ket{\psi_{\mathbf{x}^{\bar{A}}}}=\Pr(\mathbf{b}^{A}|\mathbf{x}^{A}).$
(91)

This completes the proof. ∎

We remark that as mentioned in the main text, one can set
$V_{\bar{A}}=\mathbb{I}_{\bar{A}}$ here, which does not change the relation in
Eq. 89, and thus the final conclusion. This is also demonstrated by Fig. 9.
Moreover, the used eigenbasis
$\ket{\phi_{\mathbf{x}}}:=\ket{\psi_{\mathbf{x}^{A}}}\otimes\ket{\psi_{\mathbf{x}^{\bar{A}}}}$
can also be simplified. To get the result in Section D.1, one just need to
find a basis on $\bar{A}$ of t replicas which is the eigenbasis of
$S_{t}^{\bar{A}}$. As mentioned in the main text, one thus can choose
$\ket{\phi_{\mathbf{x}}}=\ket{\psi_{\mathbf{x}^{A}}}\otimes\bigotimes_{i\in\bar{A}}\ket{\psi_{\mathbf{x}^{i}}}$,
which can be realized by
$\mathcal{R}_{A}\otimes\bigotimes_{i\in\bar{A}}\mathcal{R}_{(i)}$. In summary,
2 holds if one takes $V_{\bar{A}}=\mathbb{I}_{\bar{A}}$, and
$\ket{\phi_{\mathbf{x}}}=\ket{\psi_{\mathbf{x}^{A}}}\otimes\bigotimes_{i\in\bar{A}}\ket{\psi_{\mathbf{x}^{i}}}$,
with the function $f$ updated to

$\displaystyle f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})$
$\displaystyle\rightarrow\bigg{(}\bra{\psi_{\mathbf{x}^{A}}}\otimes\bigotimes_{i\in\bar{A}}\bra{\psi_{\mathbf{x}^{i}}}\bigg{)}\bigg{(}S_{t}^{A}\otimes\bigotimes_{i\in\bar{A}}S_{t}^{i}\bigg{)}\bigg{(}\ket{\psi_{\mathbf{x}^{A}}}\otimes\bigotimes_{i\in\bar{A}}\ket{\psi_{\mathbf{x}^{i}}}\bigg{)}$
$\displaystyle=f(\mathbf{x}^{A})\prod_{i\in\bar{A}}f(\mathbf{x}^{i}).$ (92)

We also remark that the results of the estimator and the variance proved in
the following Section D.2 also hold, after such changes.

###  D.2 Proof of Theorem 2

The proof of Theorem 2 is separated into two parts. Using 2, the first part in
Section D.2.1 shows that $\widehat{\tr_{\bar{A}}(\rho^{t})}$ defined in Eq. 17
reproduce $\tr_{\bar{A}}(\rho^{t})$ exactly in expectation (over the
randomness in the choice of unitary $V=V_{A}\otimes V_{\bar{A}}$, measurement
outcome $\mathbf{x}=\\{\mathbf{x}^{A},\mathbf{x}^{\bar{A}}\\}$, and the
mapping result $\mathbf{b}^{A}$ from $\mathbf{x}^{A}$). The second part in
Section D.2.2 proves the bound of $\mathrm{Var}(\widehat{o_{t}})$ in Eq. 18.

####  D.2.1 Unbiasedness of the estimator $\widehat{\tr_{\bar{A}}(\rho^{t})}$

The expectation of the estimator $\widehat{\tr_{\bar{A}}(\rho^{t})}$ is given
as

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}\widehat{\tr_{\bar{A}}(\rho^{t})}$
$\displaystyle=\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}{\mathrm{Re}}\left(f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})\right)\mathcal{M}_{A}^{-1}\left(V_{A}^{{\dagger}}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right)$
(93)
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x},\mathbf{b}^{A}}\Pr(\mathbf{x},\mathbf{b}^{A}|V){\mathrm{Re}}\left(f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})\right)\mathcal{M}_{A}^{-1}\left(V_{A}^{{\dagger}}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right)$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}^{A}}\mathcal{M}_{A}^{-1}\left(V_{A}^{{\dagger}}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right)\left[\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}^{A}|\mathbf{x}^{A}){\mathrm{Re}}\left(f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})\right)\right],$

where we used that
$\Pr(\mathbf{x},\mathbf{b}^{A}|V)=\Pr(\mathbf{x}|V)\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})$.
Since $\Pr(\mathbf{x}|V)$ and $\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})$ are real
numbers, the last term in the parentheses can be further written as

$\displaystyle{\mathrm{Re}}\left[\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})\right]={\mathrm{Re}}\left[\Pr(\mathbf{b}^{A}|V_{A})\right]=\bra{\mathbf{b}^{A}}V_{A}\left(\tr_{\bar{A}}(\rho^{t})\right)\ket{\mathbf{b}^{A}},$
(94)

where the first equality follows from 2. By plugging this relation into Eq.
93, we get

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}\widehat{\tr_{\bar{A}}(\rho^{t})}$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}^{A}}\mathcal{M}_{A}^{-1}\left(V_{A}^{{\dagger}}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right)\bra{\mathbf{b}^{A}}V_{A}\left(\tr_{\bar{A}}(\rho^{t})\right)\ket{\mathbf{b}^{A}}$
(95)
$\displaystyle=\mathcal{M}_{A}^{-1}\left[\mathbb{E}_{V_{A}\sim\mathcal{E}_{A}}\sum_{\mathbf{b}^{A}}\bra{\mathbf{b}^{A}}V_{A}\left(\tr_{\bar{A}}(\rho^{t})\right)\ket{\mathbf{b}^{A}}V_{A}^{{\dagger}}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right]$
$\displaystyle=\mathcal{M}_{A}^{-1}\left[\mathcal{M}_{A}\left(\tr_{\bar{A}}(\rho^{t})\right)\right]=\tr_{\bar{A}}(\rho^{t}).$

In conclusion, $\widehat{\tr_{\bar{A}}(\rho^{t})}$ is an unbiased estimator of
$\tr_{\bar{A}}(\rho^{t})$.

####  D.2.2 Proof of Eq. 18

The variance of $\widehat{o_{t}}$ reads

$\mathrm{Var}(\widehat{o_{t}})=\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}[\widehat{o_{t}}^{2}]-\left[\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}\widehat{o_{t}}\right]^{2}=\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}[\widehat{o_{t}}^{2}]-\left[\tr(O_{A}\rho^{t})\right]^{2}.$
(96)

By the self-adjoint property of $\mathcal{M}_{A}^{-1}$, the estimator
$\widehat{o_{t}}$ can be rewritten as

$\displaystyle\widehat{o_{t}}=\tr\\!\left[O_{A}\widehat{\tr_{\bar{A}}(\rho^{t})}\right]$
$\displaystyle={\mathrm{Re}}\left(f(\mathbf{x}^{A},\mathbf{x}^{\bar{A}})\right)\tr\left[O_{A}\mathcal{M}_{A}^{-1}\\!\left(V_{A}^{\dagger}\outerproduct{\mathbf{b}^{A}}{\mathbf{b}^{A}}V_{A}\right)\right]$
$\displaystyle={\mathrm{Re}}\left(\bra{\phi_{\mathbf{x}}}S_{t}\ket{\phi_{\mathbf{x}}}\right)\bra{\mathbf{b}^{A}}V_{A}\mathcal{M}_{A}^{-1}(O_{A})V_{A}^{\dagger}\ket{\mathbf{b}^{A}}.$
(97)

It follows that

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}\left[\widehat{o_{t}}^{2}\right]$
(98)
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x},\mathbf{b}^{A}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})\,\widehat{o_{t}}^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x},\mathbf{b}^{A}}\bra{\phi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\phi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}^{A}}}Q_{\mathbf{b}^{A}}\ket{\psi_{\mathbf{x}^{A}}}\left[{\mathrm{Re}}\left(\bra{\phi_{\mathbf{x}}}S_{t}\ket{\phi_{\mathbf{x}}}\right)\right]^{2}\bra{\mathbf{b}^{A}}V_{A}\mathcal{M}_{A}^{-1}(O_{A})V_{A}^{\dagger}\ket{\mathbf{b}^{A}}^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}^{A}}\bra{\mathbf{b}^{A}}V_{A}\mathcal{M}_{A}^{-1}(O_{A})V_{A}^{\dagger}\ket{\mathbf{b}^{A}}^{2}\underbrace{\left(\sum_{\mathbf{x}}\bra{\phi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\phi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}^{A}}}Q_{\mathbf{b}^{A}}\ket{\psi_{\mathbf{x}^{A}}}\right)}_{(*)},$

where $(a)$ follows because
$\Pr(\mathbf{x}|V)=\bra{\phi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\phi_{\mathbf{x}}}$ and
$\Pr(\mathbf{b}^{A}|\mathbf{x}^{A})=\bra{\phi_{\mathbf{x}}}Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes
t}\ket{\phi_{\mathbf{x}}}$, and $(b)$ follows because
$\left[\mathrm{Re}(f(\mathbf{x}))\right]^{2}=\left[\mathrm{Re}(\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}})\right]^{2}\leq
1$. The term $(*)$ can be further written as

$\displaystyle(*)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\sum_{\mathbf{x}}\bra{\phi_{\mathbf{x}}}V(\rho)^{\otimes
t}\left(Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes
t}\right)\ket{\phi_{\mathbf{x}}}\stackrel{{\scriptstyle(b)}}{{=}}\tr\left[V(\rho)^{\otimes
t}\left(Q_{\mathbf{b}^{A}}\otimes\mathbb{I}_{\bar{A}}^{\otimes
t}\right)\right]=\bra{\mathbf{b}^{A}}V_{A}(\tr_{\bar{A}}\rho)\ket{\mathbf{b}^{A}},$
(99)

where $(a)$ holds because $Q_{\mathbf{b}^{A}}$ commutes with
$\outerproduct{\phi_{\mathbf{x}}}{\phi_{\mathbf{x}}}$, $(b)$ holds because
$\\{\ket{\phi_{\mathbf{x}}}\\}_{\mathbf{x}}$ forms an orthonormal basis on
$\mathcal{H}_{d}^{\otimes t}$.

Equations (98) and (99) together imply that

$\mathbb{E}_{\\{V,\mathbf{x},\mathbf{b}^{A}\\}}\left[\widehat{o_{t}}^{2}\right]\leq\mathbb{E}_{V_{A}\sim\mathcal{E}_{A}}\sum_{\mathbf{b}^{A}}\bra{\mathbf{b}^{A}}V_{A}\mathcal{M}_{A}^{-1}(O_{A})V_{A}^{\dagger}\ket{\mathbf{b}^{A}}^{2}\bra{\mathbf{b}^{A}}V_{A}(\tr_{\bar{A}}\rho)\ket{\mathbf{b}^{A}}=\mathbb{E}\left[\tr(O_{A}\widehat{\tr_{\bar{A}}\rho})^{2}\right],$
(100)

where $\widehat{\tr_{\bar{A}}\rho}$ is the original single-copy shadow
snapshot [19] . Hence, the expectation value of $\widehat{o_{t}}^{2}$ is not
larger than that of $\left[\tr(O_{A}\widehat{\tr_{\bar{A}}\rho})\right]^{2}$
in the original shadow protocol on a single-copy of the state
$\tr_{\bar{A}}\rho$ using the unitary ensemble $\mathcal{E}_{A}$, no matter
what value $t$ takes.

By substituting Eq. 100 into Eq. 96, we have

$\displaystyle\mathrm{Var}(\widehat{o_{t}})$
$\displaystyle\leq\mathbb{E}\left[\tr(O_{A}\widehat{\tr_{\bar{A}}\rho})^{2}\right]-\left[\tr(O_{A}\rho^{t})\right]^{2}$
(101)
$\displaystyle=\mathbb{V}^{\mathcal{E}_{A}}(O_{A},\tr_{\bar{A}}\rho)+\left[\tr(O_{A}\tr_{\bar{A}}\rho)\right]^{2}-\left[\tr(O_{A}\rho^{t})\right]^{2}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{V}^{\mathcal{E}_{A}}(O_{A,0},\tr_{\bar{A}}\rho)+\left[\tr(O_{A}\rho)\right]^{2}-\left[\tr(O_{A}\rho^{t})\right]^{2}$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\left\|O_{A,0}\right\|_{\mathrm{sh},\mathcal{E}_{A}}^{2}+\|O\|_{\infty}^{2}-\left[\tr(O_{A}\rho^{t})\right]^{2},$

which confirms Eq. 18. Here,
$\mathbb{V}^{\mathcal{E}_{A}}(O_{A},\tr_{\bar{A}}\rho)$ is the variance of
$\tr(O_{A}\widehat{\tr_{\bar{A}}\rho})$ associated with the original shadow
protocol on a single-copy of $\tr_{\bar{A}}\rho$ [19], $(a)$ holds because
shifting the operator $O_{A}$ to its traceless part $O_{A,0}$ does not change
the variance of $\tr(O_{A}\widehat{\tr_{\bar{A}}\rho})$, $(b)$ holds because
$\mathbb{V}^{\mathcal{E}_{A}}(O_{A,0},\tr_{\bar{A}}\rho)\leq\left\|O_{A,0}\right\|_{\mathrm{sh},\mathcal{E}_{A}}^{2}$
[19] and $\left|\tr(O_{A}\rho)\right|\leq\|O\|_{\infty}$.

###  D.3 The estimator of moment $P_{t}$ and its variance

As mentioned in the main text, if we take $A=\emptyset$, Eq. 17 becomes

$\displaystyle\widehat{\tr(\rho^{t})}=\mathrm{Re}(f(\mathbf{x})).$ (102)

Note that such estimator can also be obtained from Eq. 4 where a global
entangling $\mathcal{R}$ is used. However, as shown in the main text, here one
can use local circuit $\bigotimes_{i=1}^{n}\mathcal{R}_{(i)}$ which entangles
the $t$ replicas in a qubit-wise manner. In this way, the eigenbasis of
$S_{t}$ is in the product form
$\ket{\phi_{\mathbf{x}}}=\bigotimes_{i=1}^{n}\ket{\psi_{\mathbf{x}^{i}}}$. As
a result, the function in Eq. 102 becomes
$f(\mathbf{x})=\prod_{i=1}^{n}\bra{\psi_{\mathbf{x}^{i}}}S_{t}^{i}\ket{\psi_{\mathbf{x}^{i}}}=\prod_{i=1}^{n}f(\mathbf{x}^{i})$,
which proves Eq. 19.

For the variance of $\widehat{\tr(\rho^{t})}$, it follows directly the general
proof of Eq. 18 in Section D.2.2. For the case $A=\emptyset$, one has
$O=\mathbb{I}$ and $O_{A,0}=0$. Inserting them in the forth line of Eq. 101,
one directly gets $\mathrm{Var}(\widehat{P_{t}})\leq
1-\left[\tr(\rho^{t})\right]^{2}$.

In Fig. 10, we numerically investigate the error of estimating the purity
$P_{2}=\tr(\rho^{2})$ using both the OS and local-AFRS protocols. Our
numerical results indicate that the estimation error of the local-AFRS
protocol diminishes as the purity of the processed quantum state increases.
Furthermore, such error is notably lower compared with that of the OS
protocol, particularly in scenarios where the purity is high. Specifically,
when $\rho$ is a pure state, the error of the local-AFRS protocol is zero,
which is consistent with our theoretical analysis outlined above.
Additionally, as shown in Fig. 10 (b), we find that the estimation error of
the local-AFRS protocol is independent of the qubit number $n$, whereas the
error of the OS protocol scales exponentially with $n$.

Figure 10: The estimation error of $P_{2}=\tr(\rho^{2})$ with respect to the
noise parameter $p$ and the qubit number $n$ for OS and local-AFRS protocols.
In (a), the quantum state is a noisy $3$-qubit GHZ state
$\rho=(1-p)\ket{\mathrm{GHZ}}\bra{\mathrm{GHZ}}+p{\mathbb{I}}_{d}/d$, in which
$p$ represents the noise parameter. In (b), the quantum state is a noisy
$n$-qubit state with $p=0.3$. The random unitary $V$ is sampled from the local
Clifford ensemble $V\in\mathcal{E}_{\mathrm{LCl}}$. We take the repetition
number $M=50$ for the local-AFRS protocol, and take $M=100$ for the OS
protocol.

##  Appendix E New estimator $\widehat{\rho^{t}}_{*}$ and its performance

In the following 1, we show a new estimator $\widehat{\rho^{t}}_{*}$ of
$\rho^{t}$ with its performance, which is constructed from Eq. 4 of Theorem 1
in the main text by averaging the estimator $\mathbf{b}$ from
$\mathbf{x}=\\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{t}\\}$.

######  Proposition 1.

Suppose one conducts the circuit proposed in Fig. 1 for once (i.e., single-
shot, $M=1$), the unbiased estimator of $\rho^{t}$ shows

$\widehat{\rho^{t}}_{*}:=\mathrm{Re}(f(\mathbf{x}))\cdot\sum_{\mathbf{b}}\Pr(\mathbf{b}|\mathbf{x})\cdot\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right),$
(103)

which satisfies
$\mathbb{E}_{\\{V,\mathbf{x}\\}}\widehat{\rho^{t}}_{*}=\rho^{t}$. Here, the
inverse channel $\mathcal{M}^{-1}$ depends on the random unitary ensemble
$\mathcal{E}$; $\mathbf{x}$ is the actual measurement outcome from the quantum
circuit in Fig. 1; the function $f$ and the mapping probability
$\Pr(\mathbf{b}|\mathbf{x})$ are given by 1.

The unbiased estimator of some observable $O$ on $\rho^{t}$, $\tr(O\rho^{t})$,
is $\widehat{o_{t}}^{*}=\tr(O\widehat{\rho^{t}}_{*})$, and its variance shows

$\displaystyle\mathrm{Var}\left(\widehat{o_{t}}^{*}\right)\leq\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)+\left[\tr(O\rho)\right]^{2},$
(104)

with

$\displaystyle\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)$
$\displaystyle:=\frac{t-1}{t}\mathbb{V}^{\mathcal{E}}_{*}(O,\rho)+\frac{1}{t}\,\mathbb{V}^{\mathcal{E}}(O,\rho),$
$\displaystyle\mathbb{V}^{\mathcal{E}}_{*}(O,\rho)$
$\displaystyle:=\underset{V\sim\mathcal{E}}{\mathbb{E}}\left(\sum_{\mathbf{b}}\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}\bra{\mathbf{b}}V(\rho)\ket{\mathbf{b}}\right)^{2}-\left[\tr(O\rho)\right]^{2},$
$\displaystyle\mathbb{V}^{\mathcal{E}}(O,\rho)$
$\displaystyle:=\underset{V\sim\mathcal{E}}{\mathbb{E}}\left(\sum_{\mathbf{b}}\bra{\mathbf{b}}V(\rho)\ket{\mathbf{b}}\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}^{2}\right)-\left[\tr(O\rho)\right]^{2}.$
(105)

Here, the value of $\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)$ is not greater than
$\mathbb{V}^{\mathcal{E}}(O,\rho)$, and $\mathbb{V}^{\mathcal{E}}(O,\rho)$ is
the variance of $\tr(O\hat{\rho})$ associated with the original shadow
protocol on a single-copy of $\rho$ [19].

We remark that $\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)$ is a quantity associated
with _multi-shot estimation_ (see e.g., Lemma 1 in [46]), which is a variant
of the original shadow protocol [19]. Such multi-shot quantity emerges here in
our protocol, as we operate same random unitary on $t$ replicas of $\rho$ in
AFRS framework. The value of $\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)$ depends on
the traceless part $O_{0}$ of the observable $O$ [47], and is not greater than
$\mathbb{V}^{\mathcal{E}}(O,\rho)$ [46]. In addition, it has been
systematically studied for various different $\mathcal{E}$ [46, 47]. In
particular, when $\mathcal{E}$ is the local Clifford ensemble
$\mathcal{E}_{\mathrm{LCl}}$ and $O=P$ is a Pauli observable with weight $w$,
we have [47, Theorem 2]

$\displaystyle\mathbb{V}^{\mathcal{E}}_{t}(P,\rho)=\frac{3^{w}}{t}+\frac{t-1}{t}\,3^{w}\left[\tr(P\rho)\right]^{2}-\left[\tr(P\rho)\right]^{2};$
(106)

when $\mathcal{E}$ is a unitary 4-design, we have [46, Theorem 2]

$\displaystyle\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)=\frac{\mathbb{V}^{\mathcal{E}}(O,\rho)}{t}+\frac{t-1}{t}\,\mathcal{O}\\!\left(2^{-n}\Tr(O_{0}^{2})\right),$
(107)

where $n$ is the number of qubits; when $\mathcal{E}$ is the ensemble of
quantum circuits constructed by interpolating $k$ single-qubit $T$ gates into
$k+1$ random Clifford layers, we have [46, Theorem 4]

$\displaystyle\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)=\frac{\mathbb{V}^{\mathcal{E}}(O,\rho)}{t}+\frac{t-1}{t}\,\mathcal{O}\\!\left(2^{-n}\Tr(O_{0}^{2})\right)+\frac{t-1}{t}\,30\Tr(O_{0}^{2})\left(1+\mathcal{O}(2^{-n})\right)\left(\frac{3}{4}+\mathcal{O}\left(2^{-n}\right)\right)^{k}.$
(108)

###### Proof of 1.

The unbiasedness of the estimator $\widehat{\rho^{t}}_{*}$ can be shown
similarly as that of $\widehat{\rho^{t}}$ in Section A.2.1. The difference is
just that the summation on the index $\mathbf{b}$ already appears in the
estimator here.

$\displaystyle\mathbb{E}_{\\{V,\mathbf{x}\\}}\widehat{\rho^{t}}_{*}$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\left[\mathrm{Re}(f(\mathbf{x}))\cdot\sum_{\mathbf{b}}\Pr(\mathbf{b}|\mathbf{x})\cdot\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)\right]$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right){\mathrm{Re}}\left[\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\Pr(\mathbf{b}|\mathbf{x})\,f(\mathbf{x})\right]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\mathcal{M}^{-1}\left(V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right)\bra{\bf{b}}V(\rho^{t})\ket{\bf{b}}$
$\displaystyle=\mathcal{M}^{-1}\left[\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}}\bra{\bf{b}}V(\rho^{t})\ket{\bf{b}}V^{{\dagger}}\outerproduct{\mathbf{b}}{\mathbf{b}}V\right]\stackrel{{\scriptstyle(b)}}{{=}}\mathcal{M}^{-1}\left[\mathcal{M}(\rho^{t})\right]=\rho^{t},$
(109)

where $(a)$ follows from 1, and $(b)$ follows from the definition of the
channel $\mathcal{M}$.

In the following, we prove the bound of
$\mathrm{Var}\big{(}\widehat{o_{t}}^{*}\big{)}$ in Eq. 104. Since the
estimator

$\displaystyle\widehat{o_{t}}^{*}=\tr\left(O\widehat{\rho^{t}}_{*}\right)=\mathrm{Re}(f(\mathbf{x}))\sum_{\mathbf{b}}\Pr(\mathbf{b}|\mathbf{x})\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}$
(110)

is a real number, its variance reads

$\displaystyle\mathrm{Var}(\widehat{o_{t}}^{*})$
$\displaystyle=\mathbb{E}_{\\{V,\mathbf{x}\\}}\left[(\widehat{o_{t}}^{*})^{2}\right]-\left(\mathbb{E}_{\\{V,\mathbf{x}\\}}\widehat{o_{t}}^{*}\right)^{2}$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\,(\widehat{o_{t}}^{*})^{2}-\left[\tr(O\rho^{t})\right]^{2}$
$\displaystyle\leq\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x}}\Pr(\mathbf{x}|V)\,\mathrm{Re}(f(\mathbf{x}))^{2}\left(\sum_{\mathbf{b}}\Pr(\mathbf{b}|\mathbf{x})\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}\right)^{2}$
$\displaystyle\leq\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}\left(\sum_{\mathbf{b}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}}\ket{\psi_{\mathbf{x}}}\bra{\mathbf{b}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}}\right)^{2}$
$\displaystyle=\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}_{1},\mathbf{b}_{2}}\bra{\mathbf{b}_{1}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}_{1}}\bra{\mathbf{b}_{2}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}_{2}}$
$\displaystyle\qquad\qquad\qquad\times\underbrace{\left(\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}_{1}}\ket{\psi_{\mathbf{x}}}\bra{\psi_{\mathbf{x}}}Q_{\mathbf{b}_{2}}\ket{\psi_{\mathbf{x}}}\right)}_{(*)},$
(111)

where the first inequality is due to ignoring the second term
$\left[\tr(O\rho^{t})\right]^{2}$, and the second inequality holds because
$\left[\mathrm{Re}(f(\mathbf{x}))\right]^{2}=\left[\mathrm{Re}(\bra{\psi_{\mathbf{x}}}S_{t}\ket{\psi_{\mathbf{x}}})\right]^{2}\leq
1$. Since $Q_{\mathbf{b}}$ and
$\outerproduct{\psi_{\mathbf{x}}}{\psi_{\mathbf{x}}}$ commute, the term $(*)$
can be written as

$\displaystyle(*)$
$\displaystyle=\sum_{\mathbf{x}}\bra{\psi_{\mathbf{x}}}V(\rho)^{\otimes
t}Q_{\mathbf{b}_{1}}Q_{\mathbf{b}_{2}}\ket{\psi_{\mathbf{x}}}\stackrel{{\scriptstyle(a)}}{{=}}\tr\\!\left[V(\rho)^{\otimes
t}Q_{\mathbf{b}_{1}}Q_{\mathbf{b}_{2}}\right]$
$\displaystyle=\frac{1}{t^{2}}\tr\\!\left[V(\rho)^{\otimes
t}\sum_{i,j=1}^{t}\left(\outerproduct{\mathbf{b}_{1}}{\mathbf{b}_{1}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)}\right)\\!\left(\outerproduct{\mathbf{b}_{2}}{\mathbf{b}_{2}}_{j}\otimes\mathbb{I}_{d}^{\otimes(t-1)}\right)\right]$
$\displaystyle=\frac{1}{t^{2}}\tr\\!\left[V(\rho)^{\otimes t}\left(\sum_{i\neq
j}\outerproduct{\mathbf{b}_{1}}{\mathbf{b}_{1}}_{i}\otimes\outerproduct{\mathbf{b}_{2}}{\mathbf{b}_{2}}_{j}\otimes\mathbb{I}_{d}^{\otimes(t-2)}+\delta_{\mathbf{b}_{1},\mathbf{b}_{2}}\sum_{i=1}^{t}\outerproduct{\mathbf{b}_{1}}{\mathbf{b}_{1}}_{i}\otimes\mathbb{I}_{d}^{\otimes(t-1)}\right)\right]$
$\displaystyle=\frac{t-1}{t}\bra{\mathbf{b}_{1}}V(\rho)\ket{\mathbf{b}_{1}}\bra{\mathbf{b}_{2}}V(\rho)\ket{\mathbf{b}_{2}}+\frac{1}{t}\,\delta_{\mathbf{b}_{1},\mathbf{b}_{2}}\bra{\mathbf{b}_{1}}V(\rho)\ket{\mathbf{b}_{1}},$
(112)

where $(a)$ holds because $\\{\ket{\psi_{\mathbf{x}}}\\}_{\mathbf{x}}$ forms
an orthonormal basis on $\mathcal{H}_{d}^{\otimes t}$. By inserting Appendix E
into Appendix E, we obtain

$\displaystyle\mathrm{Var}(\widehat{o_{t}}^{*})$
$\displaystyle\leq\frac{t-1}{t}\,\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}_{1},\mathbf{b}_{2}}\bra{\mathbf{b}_{1}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}_{1}}\bra{\mathbf{b}_{2}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}_{2}}\bra{\mathbf{b}_{1}}V(\rho)\ket{\mathbf{b}_{1}}\bra{\mathbf{b}_{2}}V(\rho)\ket{\mathbf{b}_{2}}$
$\displaystyle\quad\
+\frac{1}{t}\,\mathbb{E}_{V\sim\mathcal{E}}\sum_{\mathbf{b}_{1},\mathbf{b}_{2}}\bra{\mathbf{b}_{1}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}_{1}}\bra{\mathbf{b}_{2}}V\mathcal{M}^{-1}(O)V^{{\dagger}}\ket{\mathbf{b}_{2}}\delta_{\mathbf{b}_{1},\mathbf{b}_{2}}\bra{\mathbf{b}_{1}}V(\rho)\ket{\mathbf{b}_{1}}$
$\displaystyle=\frac{t-1}{t}\left[\mathbb{V}^{\mathcal{E}}_{*}(O,\rho)+\tr(O\rho)^{2}\right]+\frac{1}{t}\left[\mathbb{V}^{\mathcal{E}}(O,\rho)+\tr(O\rho)^{2}\right]$
$\displaystyle=\mathbb{V}^{\mathcal{E}}_{t}(O,\rho)+\left[\tr(O\rho)\right]^{2},$
(113)

which confirms Eq. 104 and completes the proof. ∎

We remark such multi-shot-like construction can also be applied to local-AFRS.
Based on Theorem 2, one has the new estimator and its variance as follows.

######  Proposition 2.

Suppose one conducts the circuit in Fig. 4 (a) for once (single-shot, $M=1$),
the unbiased estimator of $\rho^{t}$ reduced on subsystem $A$ shows