# Retrieving non-linear features from noisy quantum states

Benchi Zhao  Thrust of Artificial Intelligence, Information Hub, Hong Kong
University of Science and Technology (Guangzhou), Nansha, China  Mingrui Jing
Thrust of Artificial Intelligence, Information Hub, Hong Kong University of
Science and Technology (Guangzhou), Nansha, China  Lei Zhang  Thrust of
Artificial Intelligence, Information Hub, Hong Kong University of Science and
Technology (Guangzhou), Nansha, China  Xuanqiang Zhao  QICI Quantum
Information and Computation Initiative, Department of Computer Science, The
University of Hong Kong, Pokfulam Road, Hong Kong, China  Kun Wang  Institute
for Quantum Computing, Baidu Research, Beijing 100193, China  Xin Wang
<EMAIL_ADDRESS>Thrust of Artificial Intelligence, Information
Hub, Hong Kong University of Science and Technology (Guangzhou), Nansha, China

###### Abstract

Accurately estimating high-order moments of quantum states is an elementary
precondition for many crucial tasks in quantum computing, such as entanglement
spectroscopy, entropy estimation, spectrum estimation and predicting non-
linear features from quantum states. But in reality, inevitable quantum noise
prevents us from accessing the desired value. In this paper, we address this
issue by systematically analyzing the feasibility and efficiency of extracting
high-order moments from noisy states. We first show that there exists a
quantum protocol capable of accomplishing this task if and only if the
underlying noise channel is invertible. We then establish a method for
deriving protocols that attain optimal sample complexity using quantum
operations and classical post-processing only. Our protocols, in contrast to
conventional ones, incur lower overheads and avoid sampling different quantum
operations due to a novel technique called observable shift, making the
protocols strong candidates for practical usage on current quantum devices.
The proposed method also indicates the power of entangled protocols in
retrieving high-order information, whereas in the existing methods,
entanglement does not help. Our work contributes to a deeper understanding of
how quantum noise could affect high-order information extraction and provides
guidance on how to tackle it.

##  I Introduction

Quantum computing has emerged as a rapidly evolving field with the potential
to revolutionize the way we process and analyze information. Such an advanced
computational paradigm stores and manipulates information in a quantum state,
which forms an elaborate representation of a many-body quantum system Bennett
and DiVincenzo (2000). One critical task for this purpose is to estimate the
$k$-th moment of a quantum state’s density matrix $\rho$, which is often
denoted as $\operatorname{Tr}[\rho^{k}],k\in\mathbb{Z}^{+}$. For example, the
second moment of $\rho$ is commonly known as the purity of $\rho$. Accurately
computing $\operatorname{Tr}[\rho^{k}]$ provides an elementary precondition
for extracting spectral information of the quantum state Nielsen and Chuang
(2010), which is crucial in supporting the evaluation of non-linear functions
in quantum algorithms Childs  _et al._ (2017); Chen  _et al._ (2021), applying
to entanglement spectroscopy by determining measures of entanglement, e.g.,
Rényi entropy and von Neumann entropy Johri  _et al._ (2017); Chung  _et al._
(2014), and characterizing non-linear features of complex quantum systems in
materials Pollmann  _et al._ (2010); Vidal  _et al._ (2003); Subasi  _et al._
(2019); Li and Haldane (2008). In particular, as a core-induced development,
understanding and controlling quantum entanglement inspire various quantum
information breakthroughs including fundamental entanglement theories, quantum
cryptography, teleportation and discrimination Horodecki  _et al._ (2009); Yin
_et al._ (2020); Pirandola  _et al._ (2015); Hayashi  _et al._ (2006).

Figure 1: The general framework of recovering the high-order quantum
information $\operatorname{Tr}[\rho^{k}]$ given copies of noisy resource
${\cal N}(\rho)$ based on our derived protocol, i.e., a quantum channel ${\cal
C}$ and measurement-based post-processing. The information can be further
employed in various applications in practical quantum computing.

Numerous methods have been proposed for efficiently estimating quantum state
spectra on a quantum computer, including the deterministic quantum schemes
processing intrinsic information of the state Subramanian and Hsieh (2021);
Wang  _et al._ (2022) and the variational quantum circuit learning for
approximating non-linear quantum information functions Mitarai  _et al._
(2018); Tan and Volkoff (2021). Meanwhile, a direct estimation method of
$\operatorname{Tr}[\rho^{k}]$ through the Newton-Girard method and Hadamard-
Test Johri  _et al._ (2017) has been proposed in Subasi  _et al._ (2019), and
then it was further improved by Yirka and Subasi (2021).

However, quantum systems are inherently prone to the effects of noise, which
can arise due to a variety of factors, such as imperfect state preparation,
coupling to the environment, and imprecise control of quantum operations Clerk
_et al._ (2010). In definition, quantum noise can be described in a language
of quantum operation denoted as ${\cal N}$. Such an operation can inevitably
pose a significant challenge to the reliable estimation of
$\operatorname{Tr}[\rho^{k}]$ from corrupted copies of quantum state ${\cal
N}(\rho)$.

Previous works concentrated on the first order situation by applying the
inverse operation ${\cal N}^{-1}$ Temme  _et al._ (2017); Jiang  _et al._
(2021); Endo  _et al._ (2018) to each copy of the noisy state, such that
${\cal N}^{-1}\circ{\cal N}={\rm id}$, where id means identity map. Such
inverse operation might not be physically implementable, which requires the
usage of the quasi-probability decomposition (QPD) and sampling techniques,
decomposing ${\cal N}^{-1}$ into a linear combination of quantum channels
${\cal N}^{-1}=\sum_{i}c_{i}{\cal C}_{i}$. Then, the value
$\operatorname{Tr}[O\rho]$ can be estimated in a statistical manner, and the
total required sampling times are square proportional to sampling overhead
$g=\sum_{i}|c_{i}|$ Hoeffding (1994). Nevertheless, the situations for
estimating $\operatorname{Tr}[\rho^{k}]$ with $k>1$ stay unambiguous apart
from handling individual state noise. In this paper, we are going to retrieve
the $k$-th moment from noisy states, which is illustrated in Fig. 1. To
systematically analyze the feasibility and efficiency of extracting high-order
moment information from noisy states, as shown in Fig. 1, The following two
questions are addressed:

  1. 1.

Under what conditions can we retrieve the high-order moments from noisy
quantum states?

  2. 2.

For such conditions, what is the quantum protocol that achieves the optimal
sampling complexity?

These two questions address the existence and efficiency of quantum protocols
for retrieving high-order moment information and essential properties from
noisy states, which help us to access accurate non-linear feature estimations.

In the present study, we aim to address both of these questions. For the first
question, we establish a necessary and sufficient condition for the retrieval
of high-order moments from noisy states, which states that a quantum protocol
can achieve this goal if and only if the noisy channel is invertible.
Regarding the second question, we propose a quantum protocol that can attain
optimal sampling complexity using quantum operations and classical post-
processing only. In contrast to the conventional sampling techniques, our
protocol only employs one quantum operation due to avoiding quasi-probability
decomposition and developing a novel technique called observable shift.

We also demonstrate the advantages of our method over existing QPD methods
Temme  _et al._ (2017); Endo  _et al._ (2018) with step-by-step protocols for
some types of noise of common interest. Our protocols incur lower sampling
overheads and have simple workflows, serving as strong candidates for
practical usage on current quantum devices. The proposed method also indicates
the power of entanglement in retrieving high-order information, whereas in the
existing methods, entangled protocols do not help Jiang  _et al._ (2021);
Regula  _et al._ (2021). In the end, numerical experiments are performed to
demonstrate the effectiveness of our protocol with depolarizing noise applied
on the ground state of the Fermi-Hubbard model. Our sampling results
illustrate a more accurate estimation on $\operatorname{Tr}[\rho^{2}]$
compared with no protocol applied.

##  II Moment recoverability

In this section, we are going to address the first question proposed in the
introduction. We discover a necessary and sufficient condition for the
existence of a high-order moment extraction protocol as shown in Theorem 1.

###### Theorem 1

(Necessary and sufficient condition for existence of protocol) Given a noisy
channel ${\cal N}$, there exists a quantum protocol to extract the $k$-th
moment $\operatorname{Tr}[\rho^{k}]$ for any state $\rho$ if and only if the
noisy channel ${\cal N}$ is invertible.

Intuitively, we can understand Theorem 1 from the following aspects.
Estimating high-order moment demands complete information about quantum
channels. If a noise channel ${\cal N}$ is invertible, it means information
stored in quantum states is deformed, which can be carefully re-deformed back
to the original information with extra resources of noisy states and sampling
techniques. However, when the loss of information is unattainable, i.e., the
noise is non-invertible. Part of the information stored in the quantum state
is destroyed completely, leading to an infeasible estimation problem even with
extra quantum resources. An illustration of the theorem is shown in Fig. 2.

Figure 2: Illustration of Theorem 1. Suppose a state $\rho$ is corrupted by an
invertible channel ${\cal N}$, and $H$ is the moment observable, such that
$\operatorname{Tr}[H\rho^{\otimes k}]=\operatorname{Tr}[\rho^{k}]$. The state
information is deformed but can be retrieved via applying ${\cal D}$ and post-
processing (Top). However, for non-invertible ${\cal N}$, the high-order
moment is completely destroyed and cannot be retrieved (Bottom).

In the following, we will present a sketch of proof for our main theorem.
Starting with the definition of a quantum protocol which is usually described
as a sequence of realizable quantum operations and post-processing steps used
to perform a specific task in the domain of quantum information processing.
Mathematically, we say there exists a quantum protocol to retrieve the $k$-th
moment from copies of a noisy state ${\cal N}(\rho)$ if there exists an
operation ${\cal D}$ such that

$\operatorname{Tr}[H{\cal D}\circ{\cal N}^{\otimes k}(\rho^{\otimes
k})]=\operatorname{Tr}[H\rho^{\otimes k}],$ (1)

where $H$ is what we call the moment observable, as the usage of it is the
core of extracting the high-order moment from quantum states, i.e.,
$\operatorname{Tr}[H\rho^{\otimes k}]=\operatorname{Tr}[\rho^{k}]$. For
example, in estimating the purity of single-qubit states, the moment
observable $H$ is just a SWAP operator correlating two qubits. It is proved in
the appendix that for any order $k$, there exists such a moment observable $H$
to extract the $k$-th moment information. The inspiration of our proof comes
from the QPD method used to simulate Hermitian-preserving maps on quantum
devices, which has enjoyed great success in a variety of tasks, such as error
mitigation Temme  _et al._ (2017); Endo  _et al._ (2018), and entanglement
detection Peres (1996); Horodecki  _et al._ (2009). We extend our allowed
operation ${\cal D}$ to the field covering the Hermitian-preserving maps.

If the noisy channel ${\cal N}$ is invertible, then there exists the inverse
operation of the noisy channel ${\cal N}^{-1}$, which is generally a
Hermitian-preserving map, and $({\cal N}^{-1})^{\otimes k}$ stands for a
feasible solution to the high-order moment retriever. On the other hand, by
assuming a Hermitian-preserving map ${\cal D}$ satisfying Eq. (1) and non-
invertible ${\cal N}$. In the view of the Heisenberg picture, the adjoint of
the maps in Eq. (1) satisfies:

$\operatorname{Tr}[({\cal N}^{\otimes k})^{\dagger}\circ{\cal
D}^{\dagger}(H)\rho^{\otimes k}]=\operatorname{Tr}[H\rho^{\otimes k}].$ (2)

It has been proved in Zhao  _et al._ (2023) that given an observable $O$, a
Hermitian-preserving map ${\cal M}$ satisfies $\operatorname{Tr}[{\cal
M}(\rho)O]=\operatorname{Tr}[\rho O]$ for any state $\rho$ if and only if it
holds that ${\cal M}^{\dagger}(O)=O$. Thus, we can derive that as long as we
find a Hermitian-preserving operation ${\cal D}$ such that the condition

$({\cal N}^{\otimes k})^{\dagger}\circ{\cal D}^{\dagger}(H)=H$ (3)

is satisfied, the problem is solved. Since the effective rank of $H$ is full,
whose definition and proof are shown in appendix, then from the fact that
$\operatorname{Rank}(B)\leq\min\left(\operatorname{Rank}(A),\operatorname{Rank}(B)\right)$,
we can deduce that ${\cal N}$ is invertible, contradicting to our assumption.
This means there exists no quantum protocol for extracting high-order moments
when the noise is non-invertible. The detailed proof is given in the appendix.

##  III Observable shift method

In the previous part, we mentioned that applying the inverse operation of a
noisy channel ${\cal N}^{-1}$ to noisy states simultaneously to mitigate the
error is one feasible solution to retrieve high-order moments. However, this
channel inverse method requires exponentially many resources with respect to
$k$ to retrieve the $k$-th moment. Also, the implementation of inverse
operation ${\cal N}^{-1}$ is not quantum device friendly because it has to
sample and implement different quantum channels probabilistically.

In this section, we propose a new method called observable shift to retrieve
high-order moment information from noisy states, which requires only one
quantum operation with comparable sampling complexity.

###### Lemma 2

(Observable shift) Given an invertible quantum channel ${\cal N}$ and an
observable $O$, there exists a quantum channel ${\cal C}$, called retriever,
and coefficients $t,f$ such that

${\cal N}^{\dagger}\circ{\cal C}^{\dagger}(O)=\frac{1}{f}\left(O+tI\right).$
(4)

We develop this observable shift technique since the expectation of $O+tI$
regarding any quantum states can be computed as $\operatorname{Tr}[O\rho]+t$
during the measurement procedures. Moreover, if one wants to maintain the
retrievability of $\operatorname{Tr}[O\rho^{k}]$ from the noise channel ${\cal
N}$ with respect to any possible quantum states, then the only change that
could be made to the observable is to add constant identity since such a
transformation could maintain the original information of $\rho$. Therefore,
the trace value can be retrieved via measurement and post-processing. For
instance, when we estimate $\operatorname{Tr}[O\rho]$, where $O=I+X+Z$, By
skipping the identity, often called shifting the observable, the value of
$\operatorname{Tr}[O\rho]$ can still be re-derived by post-adding a value of
one to the expectation value of the shifted observable $O^{\prime}=X+Z$, i.e.,
$\operatorname{Tr}[O\rho]=1+\operatorname{Tr}[O^{\prime}\rho]$.

Besides, instead of mitigating noise states individually, our method utilizes
entanglement to retrieve the information with respect to the moment observable
$H$. Compared with the channel inverse method, the proposed observable shift
method requires fewer quantum resources, and its implementation is easier. The
proposed observable shift method leads to Proposition 3.

###### Proposition 3

Given error tolerance $\delta$, the $k$-th moment information can be retrieved
by sampling one quantum channel with complexity
$\mathcal{O}(f_{\min}^{2}({\cal N},k)/\delta^{2})$ and post-processing. The
quantity $f_{\min}({\cal N},k)$ is the sampling overhead defined as

$\displaystyle f_{\min}({\cal N},k)=\min\Big{\\{}f\mid$ $\displaystyle({\cal
N}^{\otimes k})^{\dagger}\circ{\cal
C}^{\dagger}(H)=\frac{1}{f}\left(H+tI\right),$ $\displaystyle
f\in\mathbb{R}^{+},t\in\mathbb{R},{\cal C}\in{\rm CPTP}\Big{\\}},$ (5)

where ${\cal N}$ is the noisy channel, ${\cal C}$ is quantum channel, $t$ is
the shifted distance, $H$ is the moment observable.

In our cases, we aim to find a Hermitian-preserving map ${\cal D}$ such that
Eq. (3) holds together with the allowance of observable shifting (4), i.e.,

$({\cal N}^{\otimes k})^{\dagger}\circ{\cal D}^{\dagger}(H)=H^{\prime}-tI,$
(6)

where $H^{\prime}=H+tI$ is the shifted observable, and $t$ is a real
coefficient. Note that the quantum channel ${\cal N}^{\otimes k}$ is a
completely positive and trace preserving (CPTP) map, and the adjoint of a CPTP
map is completely positive unital preserving Khatri and Wilde (2020), which
refers $({\cal N}^{\otimes k})^{\dagger}(I)=I$. Thus we have

$({\cal N}^{\otimes k})^{\dagger}({\cal D}^{\dagger}(H)+tI)=H^{\prime},$ (7)

where we can consider ${\cal D}^{\dagger}(H)+tI$ as a whole and denote it as
$\tilde{{\cal D}}(H)$. With proper coefficient $t$, the map $\tilde{{\cal D}}$
could reduce to a completely positive map ${\cal C}$. The detailed proof is
given in the appendix.

If we apply the quantum channel ${\cal C}$ to a noisy state and make
measurement over moment observable $H$, the expectation value will be

$\displaystyle\zeta$ $\displaystyle=\operatorname{Tr}[H{\cal C}\circ{\cal
N}^{\otimes k}(\rho^{\otimes k})]=\operatorname{Tr}[({\cal N}^{\otimes
k})^{\dagger}\circ{\cal C}^{\dagger}(H)\rho^{\otimes k}]$ (8)
$\displaystyle=\frac{1}{f}\operatorname{Tr}[(H+tI)\rho^{\otimes
k}]=\frac{1}{f}(\operatorname{Tr}[H\rho^{\otimes k}]+t).$ (9)

Obviously, the desired high-order moment is given by
$\operatorname{Tr}[H\rho^{\otimes k}]=f\zeta-t$. In order to obtain the target
expectation value of $\operatorname{Tr}[H\rho^{\otimes k}]$ within an error
$\delta$ with a probability no less than $1-p$, the number of total sampling
times $T$ is given by Hoeffding’s inequality Hoeffding (1994),

$T\geq f^{2}\frac{2}{\delta^{2}}\log(\frac{2}{p}).$ (10)

Usually, the success probability $1-p$ is fixed. Thus we consider it as a
constant in this paper, and the corresponding sample compleixty is ${\cal
O}(f^{2}/\delta^{2})$, which only depends on error tolerance $\delta$ and the
sampling overhead $f$. It is desirable to find a quantum retriever ${\cal C}$
and shift distance $t$ to make the sampling overhead $f$ as small as possible.
The optimal sampling overhead $f_{\min}({\cal N},k)$ of our method can be
calculated by SDP as follows:

$\displaystyle f_{\min}({\cal N},k)$ $\displaystyle=\min\quad f$ (11a) subject
to $\displaystyle\quad J_{\tilde{{\cal C}}}\geq 0$ (11b)
$\displaystyle\quad\operatorname{Tr}_{C}[J_{\tilde{{\cal C}}_{BC}}]=fI_{B}$
(11c) $\displaystyle\quad J_{{\cal
F}_{AC}}\equiv\operatorname{Tr}_{B}[(J_{{\cal N}_{AB}^{\otimes
k}}^{T_{B}}\otimes I_{C})(I_{A}\otimes J_{\tilde{{\cal C}}_{BC}})]$ (11d)
$\displaystyle\quad\operatorname{Tr}_{C}[(I_{A}\otimes H_{C}^{T})J_{{\cal
F}_{AC}}^{T}]=H_{A}+tI_{A}.$ (11e)

The $J_{\tilde{{\cal C}}}$ and $J_{{\cal N}^{\otimes k}}$ are the Choi-
Jamiołkowski matrices for the completely positive trace-scaling map
$\tilde{{\cal C}}=f{\cal C}$ and noise channel ${\cal N}^{\otimes k}$
respectively. Eq. (11b) corresponds to the condition that the map
$\tilde{{\cal C}}$ is completely positive, and Eq. (11c) guarantees that
$\tilde{{\cal C}}$ is a trace-scaling map. In Eq. (11d), $J_{\cal F}$ is the
Choi matrix of the composed map $\tilde{{\cal C}}\circ{\cal N}^{\otimes k}$.
Eq. (11e) corresponds to the constraint shown in Eq. (4).

Beyond retrieving particular non-linear features, i.e.,
$\operatorname{Tr}[\rho^{k}]$, we can also apply our method to estimate non-
linear functions. For a toy example, if we wish to estimate the function
$F(\rho)=\frac{1}{2}\operatorname{Tr}[\rho^{2}]+\frac{1}{3}\operatorname{Tr}[\rho^{3}]$,
we should design the moment observable first, which is supposed to be
$H=\frac{1}{2}H_{2}\otimes I+\frac{1}{3}H_{3}$, where $H_{2}$ and $H_{3}$ are
the moment observables for two and three qubits respectively. Then, the
retrieving protocol with optimal sampling overhead is given by SDP as shown in
Eq. (11).

##  IV Protocols for particular noise channels

We have introduce the observable shift method in previous part, next will
provide the analytical protocol for retrieving the second-order moment
information $\operatorname{Tr}[\rho^{2}]$ from noisy quantum states suffering
from depolarizing channel and amplitude dimpling channel, respectively.

Depolarizing channel has been extensively studied due to its simplicity and
ability to represent a wide range of physical processes that can affect
quantum states Nielsen and Chuang (2010). A quantum state undergoes a
depolarizing channel would be randomly replaced by a maximally mixed state
with a certain error rate. The single-qubit depolarizing (DE) noise ${\cal
N}_{\rm DE}^{\epsilon}$ has an exact form,

${\cal N}_{\rm DE}^{\epsilon}(\rho)=(1-\epsilon)\rho+\epsilon\frac{I}{2},$
(12)

where $\epsilon$ is the noise level, and $I$ refers to the identity operator.

Given many copies of such noisy quantum states, our method derives a protocol
for retrieving the second-order moment $\operatorname{Tr}[\rho^{2}]$ using
only one quantum channel and post-processing. Specifically, we have
Proposition 4.

###### Proposition 4

Given noisy states ${\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes 2}$, and error
tolerance $\delta$, the second order moment $\operatorname{Tr}[\rho^{2}]$ can
be estimated by $f\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\epsilon}_{\rm
DE}(\rho)^{\otimes 2}]-t$, with optimal sample complexity
$\mathcal{O}(1/(\delta^{2}(1-\epsilon)^{4})$, where
$f=\frac{1}{(1-\epsilon)^{2}}$,
$t=\frac{1-(1-\epsilon)^{2}}{2(1-\epsilon)^{2}}$. The term
$\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes
2}]$ can be estimated by implementing a quantum retriever ${\cal C}$ on noisy
states and making measurement over moment observable $H$. Moreover, there
exists an ensemble of unitary operations $\\{p_{j},U_{j}\\}_{j}$ such that the
action of the retriever ${\cal C}$ can be interpreted as,

${\cal C}(\cdot)=\sum_{j=1}p_{j}U_{j}(\cdot)U_{j}^{\dagger}$ (13)

We derive an explicit form of a mixed-unitary ensemble containing twelve fixed
unitary $U_{j}$’s given in the appendix.

As a result, the second order moment can be retrieved from depolarized states
${\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes 2}$ by applying the unitaries
$U_{j}$ randomly with equal probabilities and then performing measurements
with respect to the moment observable $H$. After repeating these steps for $T$
rounds, where $T$ is given by Eq. (10), and averaging the measurement results,
we can obtain the estimated expectation value $\zeta=\operatorname{Tr}[H{\cal
C}({\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes 2})]$. Then, the desired the
second-order moment is given by

$\operatorname{Tr}[\rho^{2}]=\frac{1}{(1-\epsilon)^{2}}\zeta-\frac{1-(1-\epsilon)^{2}}{2(1-\epsilon)^{2}}.$
(14)

When estimating $\operatorname{Tr}[\rho^{2}]$ from copies of the noisy state
${\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes 2}$, QPD-based methods incurs a
sampling overhead $\frac{(1+\epsilon/2)^{2}}{(1-\epsilon)^{2}}$. On the other
hand, our observable shift method offers a protocol with a lower sampling
overhead $\frac{1}{(1-\epsilon)^{2}}$, which is much lower than that of the
QPD-based method.

Besides, the quantum amplitude damping (AD) channel is another important model
that we are interested in, which often appears in superconducting qubits or
trapped ions. This type of noise is particularly relevant for the loss of
energy or the dissipation of excited states Breuer and Petruccione (2002),
whose action results in the transition of a qubit’s excited state to its
ground state, offering a more realistic representation of energy relaxation
processes in quantum systems. The AD channel ${\cal N}_{\rm AD}^{\varepsilon}$
is characterized by a single parameter $\varepsilon$, representing the damping
rate, which has two Kraus operators:
$A_{0}^{\varepsilon}\coloneqq|0\rangle\\!\langle
0|+\sqrt{1-\varepsilon}|1\rangle\\!\langle 1|$ and
$A_{1}^{\varepsilon}\coloneqq\sqrt{\varepsilon}|0\rangle\\!\langle 1|$, where
$\varepsilon\in[0,1]$. Similarly, given many copies of AD-produced quantum
states, the second-order information $\operatorname{Tr}[\rho^{2}]$ can be
retrieved by applying only one quantum channel and post-processing with our
protocol in Proposition 5.

###### Proposition 5

Given noisy states ${\cal N}^{\varepsilon}_{\rm AD}(\rho)^{\otimes 2}$, and
error tolerance $\delta$, the second order moment
$\operatorname{Tr}[\rho^{2}]$ can be estimated by $f\operatorname{Tr}[H{\cal
C}\circ{\cal N}^{\varepsilon}_{\rm AD}(\rho)^{\otimes 2}]-t$, with optimal
sample complexity $\mathcal{O}(1/(\delta^{2}(1-\varepsilon)^{4})$, where
$f=\frac{1}{(1-\varepsilon)^{2}}$,
$t=-\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}$. The term
$\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\varepsilon}_{\rm
AD}(\rho)^{\otimes 2}]$ can be estimated by implementing a quantum retriever
${\cal C}$ on noisy states and making measurement. Moreover, the Choi matrix
of such the retriever ${\cal C}$ is

$\displaystyle J_{\cal C}$ $\displaystyle=|00\rangle\\!\langle
00|\otimes\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$
$\displaystyle\quad+|\Psi^{+}\rangle\\!\langle\Psi^{+}|\otimes\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$
$\displaystyle\quad+|\Psi^{-}\rangle\\!\langle\Psi^{-}|\otimes\frac{1}{2}(II-H)$
$\displaystyle\quad+|11\rangle\\!\langle 11|\otimes\frac{1}{6}(II+H),$ (15)

where
$|\Psi^{\pm}\rangle\\!\langle\Psi^{\pm}|=\frac{1}{2}(|01\rangle\pm|10\rangle)(\langle
01|\pm\langle 10|)$ are Bell states.

The above retriever ${\cal C}$ can be implemented based on the following
measurement and post-processing. Given amplitude damping noisy states $({\cal
N}^{\varepsilon}_{\rm AD}(\rho))^{\otimes 2}$, we make measurement in the
basis ${\cal
B}=\\{|00\rangle,|\Psi^{+}\rangle,|\Psi^{-}\rangle,|11\rangle\\}$. From the
Choi matrix of the retriever ${\cal C}$, which is shown in Eq. (5), we know
that based on the obtained measurement results, the quantum system collapses
to the states $\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4}$ correspondingly,
where
$\sigma_{1}=\sigma_{2}=\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$,
$\sigma_{3}=\frac{1}{2}(II-H)$ and $\sigma_{4}=\frac{1}{6}(II+H)$. Each state
corresponds to a fixed expectation value $\operatorname{Tr}[H\sigma_{i}]$,
which can be predetermined via direct matrix calculation with fixed $H$ and
known $\sigma_{i}$. The next step is to run sufficient shots of basis-${\cal
B}$ measurements to determine the probability of measuring each basis state,
denoted as $p_{i}$, respectively. The term $\operatorname{Tr}[H{\cal
C}\circ{\cal N}^{\varepsilon}_{\rm AD}(\rho)^{\otimes 2}]$ is then given by
the estimated value $\zeta=\sum_{i=1}^{4}p_{i}\operatorname{Tr}[H\sigma_{i}]$.
The desired second-order moment is obtained by

$\operatorname{Tr}[\rho^{2}]=\frac{1}{(1-\varepsilon)^{2}}\zeta+\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}.$
(16)

More details can be found in the appendix. The sampling overhead for QPD-based
methods is $\frac{(1+\varepsilon)^{2}}{(1-\varepsilon)^{2}}$, while the
overhead incurred by our method is still as low as
$\frac{1}{(1-\varepsilon)^{2}}$, saying that our method requires fewer quantum
resources.

##  V Comparison with existing protocols

To extract high-order moment information from noisy states, one
straightforward method is to apply an inverse operation of noisy channel
${\cal N}^{-1}$ on quantum states to mitigate error, and then perform
measurement over moment observable $H$, which is
$\operatorname{Tr}[H\rho^{\otimes k}]=\operatorname{Tr}[H{\cal
N}^{-1}\circ({\cal N}(\rho))^{\otimes k}]$. However, the map ${\cal N}^{-1}$
might not be a physical quantum channel Jiang  _et al._ (2021), thus we cannot
implement it directly on quantum system. Fortunately, we can simulate such
channel by quasi-probability decomposition, which decompose such non-physical
map into a linear combination of physical quantum channels, i.e., ${\cal
N}^{-1}=\sum_{i}c_{i}{\cal C}_{i}$, where $c_{i}$ are the real coefficients
and ${\cal C}_{i}$ are physical quantum channels. We need to note that $c_{i}$
can be negative. From the aspect of physical implementation, in the $t$-th
round of total $T$ times of sampling, we first sample a quantum channel ${\cal
C}^{(t)}$ from $\\{{\cal C}_{i}\\}$ with probability $\\{|c_{i}|/g\\}$, where
$g=\sum_{i}|c_{i}|$, and apply it to noisy state ${\cal C}^{(t)}\circ{\cal
N}(\rho)$. Then we take measurement and get result $o^{(t)}$. After $T$ rounds
of sampling, we attain an estimation for the expectation value
$\zeta=\frac{g}{T}\sum_{t=1}^{T}\text{\rm
sgn}(c_{i}^{(t)})o^{(t)}=\operatorname{Tr}[O\rho]$. The total sampling times
$T$ is also given by Hoeffding’s inequality as shown in Eq. (10). The optimal
sampling overhead $g_{\min}({\cal N})$ is given by Jiang  _et al._ (2021);
Zhao  _et al._ (2023); Regula  _et al._ (2021)

$\displaystyle g_{\min}({\cal N})=\min\Big{\\{}\sum_{i}|c_{i}|\mid$
$\displaystyle{\cal N}^{-1}=\sum_{i}c_{i}{\cal C}_{i},$ $\displaystyle
c_{i}\in\mathbb{R},{\cal C}_{i}\in\text{CPTP}\Big{\\}},$ (17)

which can be obtained by SDP as displayed in the appendix.

When we apply the channel inverse method to retrieve the $k$-th moment, we
should apply the inverse operation simultaneously on $k$ quantum systems, the
corresponding optimal sampling overhead if given by $g_{\min}({\cal N},k)$,
which is

$\displaystyle g_{\min}({\cal N},k)=\min\Big{\\{}\sum_{i}|c_{i}|\mid$
$\displaystyle({\cal N}^{-1})^{\otimes k}=\sum_{i}c_{i}{\cal C}_{i},$
$\displaystyle c_{i}\in\mathbb{R},{\cal C}_{i}\in\text{CPTP}\Big{\\}},$ (18)

Figure 3: The sampling overhead respective to noise level for estimating
$\operatorname{Tr}[\rho^{3}]$ from amplitude damping noise corrupted state
${\cal N}(\rho)$. The red curve refers to the overhead for channel inverse,
the green curve stands for the information recovery method, and the blue curve
is the newly proposed method.

With Eq.(3) and Eq. (V), we can make comparison of the sampling overhead
between our method and the conventional QPD channel inverse method, which
leads to Lemma 6.

###### Lemma 6

For arbitrary invertible quantum noisy channel ${\cal N}$, and moment order
$k$, we have

$f_{\min}({\cal N},k)\leq g_{\min}({\cal N},k)$ (19)

Lemma 6 implies that in the task of extracting the $k$-th moment
$\operatorname{Tr}[\rho^{k}]$ from noisy states, the proposed method requires
fewer sampling times, i.e., consumes fewer quantum resources. The detailed
proof is displayed in appendix. It has been prove by Regula et al.in Ref.
Regula  _et al._ (2021) that the sampling overhead for simulating a trace
preserving linear map is equivalent to its diamond norm. Specifically, in the
case of inverse operation ${\cal N}^{-1}$, we have $g_{\min}({\cal N})=\|{\cal
N}^{-1}\|_{\diamond}$. Note that the diamond norm is multiplicativity with
respect to tensor product Regula  _et al._ (2021); Jiang  _et al._ (2021),
i.e., $\|({\cal N}^{-1})^{\otimes k}\|_{\diamond}=\|{\cal
N}^{-1}\|_{\diamond}^{k}$. Thus, we conclude that the optimal sampling
overhead for retrieve high-moment from noisy states increases exponential
respect to the moment order $k$, which is

$\displaystyle g_{\min}({\cal N},k)$ $\displaystyle=g_{\min}({\cal N}^{\otimes
k})=\|({\cal N}^{-1})^{\otimes k}\|_{\diamond}=\|{\cal
N}^{-1}\|_{\diamond}^{k}$ $\displaystyle=g_{\min}({\cal N})^{k}.$ (20)

In order to illustrate the advantage of the proposed method over the channel
inverse method in the term of sampling overhead, we conduct a numerical
experiment to extract the third moment $\operatorname{Tr}[\rho^{3}]$ from
amplitude damping noise channel with different noise levels. The results are
shown in Fig. 3. The red and blue curves stand for the sampling overhead for
channel inverse and observable shift method, respectively.

Compared with the QPD-based inverse operation method, the proposed observable
shift method have at least two-folded advantages. First, our method can
achieve a lower sampling overhead, i.e., $f_{\min}({\cal N},k)\leq
g_{\min}({\cal N},k)$. In previous, we have shown examples where
$f_{\min}({\cal N},k)$ strictly smaller, indicating the effectiveness of the
newly proposed observable shift technique. It is also observed that the
optimal protocol given by our method is generally an entangled one. In
contrast to the uselessness of entanglement in QPD Jiang  _et al._ (2021);
Regula  _et al._ (2021), our method demonstrates the power of entanglement for
tackling noise. Second, protocols given by our method are more hardware
friendly as they only need to repeat one fixed quantum channel, whereas the
QPD-based methods have to sample from multiple quantum channels and implement
each of them.

##  VI Application to Fermi-Hubbard model

The Fermi-Hubbard model is a key focus in condensed matter physics due to its
relevance in metal-insulator transitions and high-temperature
superconductivity Cade  _et al._ (2020); Dagotto (1994). Recent studies have
shown that entanglement spectroscopy can be utilized to extract critical
exponents and phase transitions in the Fermi-Hubbard model Kokail  _et al._
(2021); Linke  _et al._ (2018); Szasz  _et al._ (2020). As the model is
characterized by a broad range of correlated electrons, it necessitates multi-
determinant and highly accurate calculations Ferreira  _et al._ (2022); Szasz
_et al._ (2020) which hence demand ingenious methods of quantum noise control.

In a physical system such as a metallic crystal with an $n_{x}\times n_{y}$
square lattice, each lattice point, known as a site, is assigned an index. The
Hubbard model Hamiltonian takes on a fermionic form in second quantization,

$\displaystyle H_{\text{Hubbard}}$ $\displaystyle=-J\sum_{\langle
i,j\rangle,\sigma}(a^{\dagger}_{i\sigma}a_{j\sigma}+a^{\dagger}_{j\sigma}a_{i\sigma})$
(21) $\displaystyle+U\sum_{i}n_{i\uparrow}n_{i\downarrow}+H_{\text{local}},$

where $a^{\dagger}_{i\sigma},a_{i\sigma}$ are fermionic creation and
annihilation operators; $n_{i\sigma}=a^{\dagger}_{i\sigma},a_{i\sigma}$, are
the number operators; the notation $\langle i,j\rangle$ associates adjacent
sites in the $n_{x}\times n_{y}$ rectangular lattice;
$\sigma\in\\{\uparrow,\downarrow\\}$ labels the spin orbital. The first term
in Eq. (21) corresponds to the hopping term, where $J$ denotes the tunnelling
amplitude. The second term involves the on-site Coulomb repulsion, represented
by $U$. The final term in the equation defines the local potential resulting
from nuclear-electron interaction, which we have chosen to be the Gaussian
form Wecker  _et al._ (2015).

$H_{\text{local}}=\sum_{j=1}\sum_{\nu=\uparrow,\downarrow}\epsilon_{j,\nu}n_{j,\nu};\quad\epsilon_{j,\nu}=-\lambda_{\nu}e^{-\frac{1}{2}(j-m_{\nu})^{2}/\sigma_{\nu}^{2}}.$
(22)

In the following, we consider a specific 3-site (6-qubit) Fermi-Hubbard
Hamiltonian with $J=2,U=3$ and $\lambda_{\uparrow,\downarrow}=3,0.1$,
$m_{\uparrow,\downarrow}=3,3$. The standard deviation $\sigma_{v}$ for both
spin-up and -down potentials are set to $1$ guaranteeing a charge-spin
symmetry around the centre site ($i=2$) of the chain system.

The ground state entanglement spectroscopy of the model identifies the
topological-ordering signatures of the system which requires high-precision
entropy estimations over each bipartite sector of the entire system. We show
that the mitigation of the quantum noise can be achieved and therefore,
enhance the determination of $\operatorname{Tr}[{\cal N}(\rho_{A})^{2}]$, via
our proposed method, which is displayed in the appendix.

Fig. 4 displays the sampling distribution with and without error mitigation.
The orange curve refers to the estimation distribution of second-order
information $\operatorname{Tr}[\rho^{2}]$ from noisy states, and the cyan
curve shows the estimation distribution with error mitigation. And the black
dash line is the exact value of $\operatorname{Tr}[\rho^{2}]$.

Figure 4: Simulation of high-order moment $\operatorname{Tr}[\rho^{2}]$
estimation. The curves are calculated from sampling. The orange curve
represents the estimation from depolarizing noised state ${\cal N}(\rho)$ with
noise level $\epsilon=0.1$. The cyan curve is the estimation with the proposed
error mitigation method. The black dashed line stands for the exact value of
$\operatorname{Tr}[\rho^{2}]$.

##  VII Conclusion and discussion

In this study, we establish that when quantum states are distorted by noises,
the original moment information can still be retrieved through post-processing
if and only if the noise is invertible. Furthermore, our proposed method,
called observable shift, outperforms QPD-based techniques from two aspects:
(1) The proposed method requires lower quantum sampling complexity than the
existing one, which implies the superiority of entangled protocols over
product protocols. This contrasts with the multiplicativity of cost observed
in QPD-based methods for quantum error mitigation. (2) The observable shift
method is easier to implement than the QPD-based method as it only involves a
single quantum operation, which makes our method more friendly to quantum
devices. Our findings have implications for the dependable estimation of non-
linear information in quantum systems and can influence various applications,
including entanglement spectroscopy and ground state property estimation.

For further work, it will be interesting to improve the scalability of the
observable shift method, which makes this approach more practical. We expect
the observable shift technique can be incorporated into more algorithms and
protocols to boost efficiency.

##  VIII Acknowledgements

Authors would like to thank Chengkai Zhu and Chenghong Zhu for their valuable
discussion.

## References

  * Bennett and DiVincenzo (2000) C. H. Bennett and D. P. DiVincenzo, Nature 404, 247 (2000), number: 6775 Publisher: Nature Publishing Group.
  * Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Childs  _et al._ (2017) A. M. Childs, R. Kothari, and R. D. Somma, SIAM Journal on Computing 46, 1920 (2017).
  * Chen  _et al._ (2021) R. Chen, Z. Song, X. Zhao, and X. Wang, Quantum Science and Technology 7, 015019 (2021).
  * Johri  _et al._ (2017) S. Johri, D. S. Steiger, and M. Troyer, Physical Review B 96, 195136 (2017).
  * Chung  _et al._ (2014) C.-M. Chung, L. Bonnes, P. Chen, and A. M. Läuchli, Physical Review B 89, 195147 (2014).
  * Pollmann  _et al._ (2010) F. Pollmann, A. M. Turner, E. Berg, and M. Oshikawa, Phys. Rev. B 81, 064439 (2010).
  * Vidal  _et al._ (2003) G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. 90, 227902 (2003).
  * Subasi  _et al._ (2019) Y. Subasi, L. Cincio, and P. J. Coles, Journal of Physics A: Mathematical and Theoretical 52, 044001 (2019).
  * Li and Haldane (2008) H. Li and F. D. M. Haldane, Phys. Rev. Lett. 101, 010504 (2008).
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Reviews of modern physics 81, 865 (2009).
  * Yin  _et al._ (2020) J. Yin, Y.-H. Li, S.-K. Liao, M. Yang, Y. Cao, L. Zhang, J.-G. Ren, W.-Q. Cai, W.-Y. Liu, S.-L. Li, _et al._ , Nature 582, 501 (2020).
  * Pirandola  _et al._ (2015) S. Pirandola, J. Eisert, C. Weedbrook, A. Furusawa, and S. L. Braunstein, Nature photonics 9, 641 (2015).
  * Hayashi  _et al._ (2006) M. Hayashi, D. Markham, M. Murao, M. Owari, and S. Virmani, Physical review letters 96, 040501 (2006).
  * Subramanian and Hsieh (2021) S. Subramanian and M.-H. Hsieh, Physical review A 104, 022428 (2021).
  * Wang  _et al._ (2022) Y. Wang, L. Zhang, Z. Yu, and X. Wang, arXiv preprint arXiv:2209.14278 (2022).
  * Mitarai  _et al._ (2018) K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Phys. Rev. A 98, 032309 (2018).
  * Tan and Volkoff (2021) K. C. Tan and T. Volkoff, Physical Review Research 3, 033251 (2021).
  * Yirka and Subasi (2021) J. Yirka and Y. Subasi, Quantum 5, 535 (2021).
  * Clerk  _et al._ (2010) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Reviews of Modern Physics 82, 1155 (2010).
  * Temme  _et al._ (2017) K. Temme, S. Bravyi, and J. M. Gambetta, Physical review letters 119, 180509 (2017).
  * Jiang  _et al._ (2021) J. Jiang, K. Wang, and X. Wang, Quantum 5, 600 (2021).
  * Endo  _et al._ (2018) S. Endo, S. C. Benjamin, and Y. Li, Physical Review X 8, 031027 (2018).
  * Hoeffding (1994) W. Hoeffding, The collected works of Wassily Hoeffding , 409 (1994).
  * Regula  _et al._ (2021) B. Regula, R. Takagi, and M. Gu, Quantum 5, 522 (2021).
  * Peres (1996) A. Peres, Physical Review Letters 77, 1413 (1996).
  * Zhao  _et al._ (2023) X. Zhao, B. Zhao, Z. Xia, and X. Wang, Quantum 7, 978 (2023).
  * Khatri and Wilde (2020) S. Khatri and M. M. Wilde, arXiv preprint arXiv:2011.04672 (2020).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, USA, 2002).
  * Cade  _et al._ (2020) C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Physical Review B 102 (2020), 10.1103/physrevb.102.235122.
  * Dagotto (1994) E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).
  * Kokail  _et al._ (2021) C. Kokail, B. Sundar, T. V. Zache, A. Elben, B. Vermersch, M. Dalmonte, R. van Bijnen, and P. Zoller, Phys. Rev. Lett. 127, 170501 (2021).
  * Linke  _et al._ (2018) N. M. Linke, S. Johri, C. Figgatt, K. A. Landsman, A. Y. Matsuura, and C. Monroe, Physical Review A 98, 052334 (2018).
  * Szasz  _et al._ (2020) A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Physical Review X 10, 021042 (2020).
  * Ferreira  _et al._ (2022) D. L. Ferreira, T. O. Maciel, R. O. Vianna, and F. Iemini, Physical Review B 105, 115145 (2022).
  * Wecker  _et al._ (2015) D. Wecker, M. B. Hastings, and M. Troyer, Physical Review A 92 (2015), 10.1103/physreva.92.042303.

Supplementary Material and Proofs for
Retrieving non-linear features from noisy quantum states

##  Appendix A Symbols and Notations

In this work, we focus on linear maps ${\cal N}$ having the same input and
output dimension. A linear map ${\cal N}$ is called Hermitian-preserving (HP)
if, for any introduced reference system $R$, the product map
$id_{R}\otimes{\cal N}$ maps any Hermitian operator to another Hermitian
operator. If additionally the map also preserves the positivity of the
operators. We say it is completely positive (CP). Besides, a linear map ${\cal
N}$ is trace scaling (TS) means $\operatorname{Tr}[{\cal
N}(X)]=a\operatorname{Tr}[X]$ for all $X\in{\cal L}_{A},a\in\mathbb{R}$. If
$a=1$, such map ${\cal N}$ is called trace-preserving (TP). We call a linear
map ${\cal N}_{A\rightarrow A^{\prime}}$ a quantum channel if it is both
completely positive and trace preserving (CPTP). More, by saying ${\cal N}$ is
unital scaling (US), we mean $\operatorname{Tr}[{\cal
N}(I)]=b\operatorname{Tr}[I]$, where $I$ refers to the identity operator, and
$b\in\mathbb{R}$. Again, if $b=1$, we say ${\cal N}$ is unital preserving
(UP). A sketch of the relationships among these classes of linear maps has
been shown in the Venn diagram 5.

Remark: Adjoint map has the following properties Khatri and Wilde (2020); Zhao
_et al._ (2023)

  * •

The adjoint of a CP map is CP;

  * •

The adjoint of a TP map is UP;

  * •

The adjoint of a UP map is TP;

  * •

The adjoint of a TS map is US;

  * •

The adjoint of a US map is TS;

Figure 5: The Venn diagram summarizing relationships among different classes
of linear maps stated in the section of notations.

##  Appendix B Proof of the existence of observable for high-order moment

###### Lemma 7

Suppose $\rho$ is an $n$-qubit state. Then for a positive integer $k$, there
exists an $nk$-qubit Hermitian matrix $H$ such that
$\operatorname{Tr}\left[H\rho^{\otimes
k}\right]=\operatorname{Tr}\left[\rho^{k}\right]$.

###### Proof.

Consider an orthonormal basis $\\{|x\rangle\\}_{x}$ for an $n$-qubit system.
The desired Hermitian matrix can be constructed from an $nk$-qubit cyclic
permutation matrix

$S_{k}=\sum_{\bm{x}\coloneqq x_{1}\cdots
x_{k}=0}^{2^{nk}-1}\bigotimes_{j=1}^{k}|x_{\pi(j)}\rangle\\!\langle x_{j}|,$
(23)

where $\pi=(1,\ldots,k)$ is a permutation function in cyclic notation. Note
that the conjugate transpose of this matrix permutes the subsystems in an
inverse order i.e.

$S_{k}^{\dagger}=\sum_{\bm{x}}\bigotimes_{j=1}^{k}|x_{j}\rangle\\!\langle
x_{\pi(j)}|=\sum_{\bm{x}}\bigotimes_{j=1}^{k}|x_{\pi^{-1}(j)}\rangle\\!\langle
x_{j}|.$ (24)

In the rest of this proof, we will prove that

$\operatorname{Tr}\left[S_{k}\rho^{\otimes
k}\right]=\operatorname{Tr}\left[S_{k}^{\dagger}\rho^{\otimes
k}\right]=\operatorname{Tr}\left[\rho^{k}\right]$ (25)

for all $\rho$, and hence the result follows by setting
$H=\frac{1}{2}\left(S_{k}+S_{k}^{\dagger}\right)$.

To start with, decomposing $\rho=\sum_{x,y}C_{xy}|x\rangle\\!\langle y|$ gives

$\displaystyle\rho^{k}$
$\displaystyle=\sum_{\bm{x},\bm{y}}\prod_{j=1}^{k}C_{x_{j}y_{j}}|x_{j}\rangle\\!\langle
y_{j}|=\sum_{\bm{x},y_{k}}\left(\prod_{j=1}^{k-1}C_{x_{j}x_{\pi(j)}}\right)c_{x_{k}y_{k}}|x_{1}\rangle\\!\langle
y_{k}|\implies\operatorname{Tr}\left[\rho^{k}\right]=\sum_{\bm{x}}\prod_{j=1}^{k}C_{x_{j}x_{\pi(j)}},$
(26) $\displaystyle\rho^{\otimes k}$
$\displaystyle=\sum_{\bm{x},\bm{y}}\bigotimes_{j=1}^{k}C_{x_{j}y_{j}}|x_{j}\rangle\\!\langle
y_{j}|.$ (27)

Apply $S_{k}$ on $\rho^{\otimes k}$ and we have

$\displaystyle S_{k}\rho^{\otimes k}$
$\displaystyle=\sum_{\bm{x},\bm{y},\bm{z}}\bigotimes_{j=1}^{k}|z_{\pi(j)}\rangle\\!\langle
z_{j}|\cdot C_{x_{j}y_{j}}|x_{j}\rangle\\!\langle
y_{j}|=\sum_{\bm{x},\bm{y}}\bigotimes_{j=1}^{k}C_{x_{j}y_{j}}|x_{\pi(j)}\rangle\\!\langle
y_{j}|$ (28) $\displaystyle\implies\operatorname{Tr}\left[S_{k}\rho^{\otimes
k}\right]$
$\displaystyle=\sum_{\bm{x}}\prod_{j=1}^{k}C_{x_{j}x_{\pi(j)}}=\operatorname{Tr}\left[\rho^{k}\right].$
(29)

Lastly, the statement $\operatorname{Tr}\left[S_{k}^{\dagger}\rho^{\otimes
k}\right]$ holds from the fact that $\rho$ is Hermitian and thus

$\operatorname{Tr}\left[S_{k}^{\dagger}\rho^{\otimes
k}\right]=\sum_{\bm{x}}\prod_{j=1}^{k}C_{x_{j}x_{\pi^{-1}(j)}}=\sum_{\bm{x}}\prod_{j=1}^{k}C_{x_{\pi^{-1}(j)}x_{j}}^{*}=\sum_{\bm{x}}\prod_{j=1}^{k}C_{x_{j}x_{\pi(j)}}^{*}=\operatorname{Tr}\left[(\rho^{*})^{k}\right]=\operatorname{Tr}\left[\rho^{k}\right].$
(30)

$\sqcap$$\sqcup$

##  Appendix C Proof of Theorem 1

###### Lemma 8

Zhao  _et al._ (2023) Given an observable $H$, an HP map ${\cal M}$ satisfies
$\operatorname{Tr}[({\cal M}(\rho)O)]=\operatorname{Tr}[\rho O]$ for any state
$\rho$ if and only if it holds that ${\cal M}^{\dagger}(O)=O$.

###### Definition 1

(Vectorization Khatri and Wilde (2020)) The vectorization of a matrix
$M=\sum_{ij}M_{ij}|i\rangle\\!\langle j|$ is defined as

$|M\rangle\coloneqq\sum_{i,j}M_{ij}|j\rangle\otimes|i\rangle.$ (31)

###### Definition 2

(Hermitian effective rank) For a $2^{n}\times 2^{n}$ Hermitian matrix $O$ in a
bipartite system $AB$, the effective rank ${\cal R}(O)$ is defined as

${\cal R}(O)\coloneqq\operatorname{Rank}(|O\rangle\\!\langle O|_{A}),$ (32)

where $|O\rangle\\!\langle O|_{A}=\operatorname{Tr}_{B}[|O\rangle\\!\langle
O|_{AB}]$.

###### Definition 3

(Matrix of channel Zhao  _et al._ (2023)) A quantum channel ${\cal N}$ can be
represented in the matrix form $M_{\cal N}$ by its Kraus operators
$\\{E_{k}\\}$. Given ${\cal N}(|i\rangle\\!\langle
j|)=\sum_{k}E_{k}|i\rangle\\!\langle j|E_{k}^{\dagger}$, we can represent it
in vector form, which is $\text{vec}({\cal N}(|i\rangle\\!\langle
j|))=\sum_{k}\bar{E_{k}}|j\rangle\otimes
E_{k}|i\rangle=(\sum_{k}\bar{E_{j}}\otimes E_{k})|j\rangle\otimes|i\rangle$,
where $\bar{E_{k}}$ is the complex conjugate of $E_{k}$. Here we can define

$M_{\cal N}\coloneqq\sum_{k}\bar{E_{k}}\otimes E_{k}$ (33)

as the matrix representation of channel ${\cal N}$. The rank of such matrix
$M_{\cal N}$ is called the channel rank of ${\cal N}$.

###### Lemma 9

(Full effective rank) Suppose
$H=\frac{1}{2}\left(S_{k}+S_{k}^{\dagger}\right)$ is an $nk$-qubit Hermitian
matrix with $k$ subsystems. Then ${\cal R}(H)=4^{n}$ for arbitrary $n$-qubit
subsystem $A$ i.e. the effective rank of $H$ is full.

###### Proof.

By Equation (23) and Equation (24), $H$ can be decomposed as

$H=\frac{1}{2}\sum_{\bm{x}\coloneqq x_{1}\cdots
x_{k}=0}^{2^{nk}-1}\bigotimes_{j=1}^{k}|x_{\pi(j)}\rangle\\!\langle
x_{j}|+\bigotimes_{j=1}^{k}|x_{\pi^{-1}(j)}\rangle\\!\langle x_{j}|$ (34)

and thus the vectorization of $H$ is

$\displaystyle|H\rangle$
$\displaystyle=\frac{1}{2}\sum_{\bm{x}}\left(\bigotimes_{j=1}^{k}|x_{j}\rangle\right)\otimes\left(\bigotimes_{j=1}^{k}|x_{\pi(j)}\rangle\right)+\left(\bigotimes_{j=1}^{k}|x_{j}\rangle\right)\otimes\left(\bigotimes_{j=1}^{k}|x_{\pi^{-1}(j)}\rangle\right)$
(35)
$\displaystyle=\frac{1}{2}\sum_{\bm{x}}\left(\bigotimes_{j=1}^{k}|x_{j}\rangle\right)\otimes\left(\bigotimes_{j=1}^{k}|x_{\pi(j)}\rangle+\bigotimes_{j=1}^{k}|x_{\pi^{-1}(j)}\rangle\right).$
(36)

Let $A$ be the $l$-th subsystem of $H$ and $B$ as the rest of subsystems. Note
that $|H\rangle\\!\langle H|$ can be decomposed as

$\displaystyle|H\rangle\\!\langle
H|=\frac{1}{2}\sum_{\bm{x},\bm{y}}\left(\bigotimes_{j=1}^{k}|x_{j}\rangle\\!\langle
y_{j}|\right)\otimes\left(\bigotimes_{j=1}^{k}|x_{\pi(j)}\rangle+\bigotimes_{j=1}^{k}|x_{\pi^{-1}(j)}\rangle\right)\left(\bigotimes_{j=1}^{k}\langle
y_{\pi(j)}|+\bigotimes_{j=1}^{k}\langle y_{\pi^{-1}(j)}|\right).$ (37)

Applying $\operatorname{Tr}_{A}$ on $|H\rangle\\!\langle H|$ gives

$\displaystyle|H\rangle\\!\langle
H|_{A}=\frac{1}{4}\sum_{\bm{x},\bm{y}}\operatorname{Tr}\left[\bigotimes_{j\neq
l}^{k}|x_{j}\rangle\\!\langle y_{j}|\right]|x_{l}\rangle\\!\langle
y_{l}|\otimes\Bigg{\\{}$
$\displaystyle\operatorname{Tr}\left[\bigotimes_{j\neq
l}^{k}|x_{\pi(j)}\rangle\\!\langle
y_{\pi(j)}|\right]|x_{\pi(l)}\rangle\\!\langle y_{\pi(l)}|+$ (38)
$\displaystyle\operatorname{Tr}\left[\bigotimes_{j\neq
l}^{k}|x_{\pi^{-1}(j)}\rangle\\!\langle
y_{\pi^{-1}(j)}|\right]|x_{\pi^{-1}(l)}\rangle\\!\langle y_{\pi^{-1}(l)}|+$
$\displaystyle\operatorname{Tr}\left[\bigotimes_{j\neq
l}^{k}|x_{\pi(j)}\rangle\\!\langle
y_{\pi^{-1}(j)}|\right]|x_{\pi(l)}\rangle\\!\langle y_{\pi^{-1}(l)}|+$
$\displaystyle\operatorname{Tr}\left[\bigotimes_{j\neq
l}^{k}|x_{\pi^{-1}(j)}\rangle\\!\langle
y_{\pi(j)}|\right]|x_{\pi^{-1}(l)}\rangle\\!\langle y_{\pi(l)}|\Bigg{\\}}$
$\displaystyle\implies|H\rangle\\!\langle
H|_{A}=\frac{1}{4}\sum_{\bm{x},\bm{y}}\left(\prod_{j\neq
l}^{k}\delta_{x_{j}y_{j}}\right)|x_{l}\rangle\\!\langle
y_{l}|\otimes\Bigg{\\{}$ $\displaystyle\left(\prod_{j\neq
l}^{k}\delta_{x_{\pi(j)}y_{\pi(j)}}\right)|x_{\pi(l)}\rangle\\!\langle
y_{\pi(l)}|+$ (39) $\displaystyle\left(\prod_{j\neq
l}^{k}\delta_{x_{\pi^{-1}(j)}y_{\pi^{-1}(j)}}\right)|x_{\pi^{-1}(l)}\rangle\\!\langle
y_{\pi^{-1}(l)}|+$ $\displaystyle\left(\prod_{j\neq
l}^{k}\delta_{x_{\pi(j)}y_{\pi^{-1}(j)}}\right)|x_{\pi(l)}\rangle\\!\langle
y_{\pi^{-1}(l)}|+$ $\displaystyle\left(\prod_{j\neq
l}^{k}\delta_{x_{\pi^{-1}(j)}y_{\pi^{(}j)}}\right)|x_{\pi^{-1}(l)}\rangle\\!\langle
y_{\pi(l)}|\Bigg{\\}}.$

Observe that $\prod_{j\neq l}^{k}\delta_{x_{j}y_{j}}=1$ if and only if
$x_{j}=y_{j}$ for all $j\neq l$, under which case $\prod_{j\neq
l}^{k}\delta_{x_{\pi(j)}y_{\pi(j)}}=\prod_{j\neq
l}^{k}\delta_{x_{\pi(j)}y_{\pi^{-1}(j)}}=\prod_{j\neq
l}^{k}\delta_{x_{\pi^{-1}(j)}y_{\pi^{(}j)}}=\delta_{x_{l},y_{l}}$. Then

$\displaystyle|H\rangle\\!\langle H|_{A}$
$\displaystyle=\frac{1}{4}\sum_{x_{l},x_{\pi(l)},x_{\pi^{-1}(l)},y_{l}}\delta_{x_{l},y_{l}}|x_{l}\rangle\\!\langle
y_{l}|\otimes\left(|x_{\pi(l)}\rangle+|x_{\pi^{-1}(l)}\rangle\right)\left(\langle
x_{\pi(l)}|+\langle x_{\pi^{-1}(l)}|\right)$ (40)
$\displaystyle=\frac{1}{4}\left[\sum_{x_{l}=0}^{2^{n}-1}|x_{l}\rangle\\!\langle
x_{l}|\right]\otimes\left[\sum_{x_{\pi(l)},x_{\pi^{-1}(l)}=0}^{2^{n}-1}\left(|x_{\pi(l)}\rangle+|x_{\pi^{-1}(l)}\rangle\right)\left(\langle
x_{\pi(l)}|+\langle x_{\pi^{-1}(l)}|\right)\right].$ (41)

Here, the first bracket gives the identity matrix of size $2^{n}$, and the
later one gives a matrix with column space of dimension $2^{n}$. The result
follows from the fact $\operatorname{Rank}(A\otimes
B)=\operatorname{Rank}(A)+\operatorname{Rank}(B)$.  $\sqcap$$\sqcup$

Theorem 1 (Necessary and sufficient condition for existence of protocol) Given
a noisy channel ${\cal N}$, there exists a quantum protocol to extract the
$k$-th moment $\operatorname{Tr}[\rho^{k}]$ for any state $\rho$ if and only
if the noisy channel ${\cal N}$ is invertible.

###### Proof.

From Lemma 7, there exists an observable $H$ such that
$\operatorname{Tr}[H\rho^{\otimes k}]=\operatorname{Tr}[\rho^{k}]$. If the
noisy channel ${\cal N}$ is invertible, then there exists the inverse
operation of noisy channel ${\cal N}^{-1}$, such that we can apply the inverse
operation on quantum system, and perform measurement with respect to moment
observable $H$. Thus, the high-order moment can be retrieved, which is

$\displaystyle\operatorname{Tr}[H\cdot({\cal N}^{-1})^{\otimes k}\circ{\cal
N}^{\otimes k}(\rho^{\otimes k})]=\operatorname{Tr}[H\rho^{\otimes
k}]=\operatorname{Tr}[\rho^{k}].$ (42)

The proof for "if" part is completed.

For the "only if" part, let’s assume there exists a quantum protocol to
extract $k$-th order information $\operatorname{Tr}[\rho^{k}]$ from noisy
state ${\cal N}(\rho)$, with non-invertible noise channel ${\cal N}$. Then,
from the definition of quantum protocol, the existence of quantum protocol to
extract information $\operatorname{Tr}[\rho^{k}]$ refers to the existence of
an HP map ${\cal D}$ such that the following equation holds

$\operatorname{Tr}[H\cdot{\cal D}\circ{\cal N}^{\otimes k}(\rho^{\otimes
k})]=\operatorname{Tr}[H\rho^{\otimes k}]=\operatorname{Tr}[\rho^{k}].$ (43)

From Lemma 8, we convert the problem into that there exists an HP map ${\cal
D}$ such that the following equation holds

$({\cal N}^{\otimes k})^{\dagger}\circ{\cal D}^{\dagger}(H)=H.$ (44)

Since ${\cal D}$ is an HP map, then the adjoint map ${\cal D}^{\dagger}$ is an
HP map as well, thus we have $D^{\dagger}(H)=Q$, where $Q$ is hermitian. Now
we have

$({\cal N}^{\otimes k})^{\dagger}(Q)=H$ (45)

Next we vectorize the Hemritian matrix $H$ and $Q$ into $|H\rangle$ and
$|Q\rangle$, the corresponding channel $({\cal N}^{\otimes k})^{\dagger}$ in
vectoriztion representation is $M_{{\cal N}^{\dagger}}^{\otimes k}$, where
$M_{{\cal N}^{\dagger}}=\sum_{k}\bar{E_{k}}^{\dagger}\otimes E_{k}^{\dagger}$,
$\\{E_{k}\\}$ is the Kraus representation of channel ${\cal N}$ operator, and
$\bar{E_{k}}$ is the complex conjugate of $E_{k}$. Thus the equations becomes

$M_{{\cal N}^{\dagger}}^{\otimes k}|Q\rangle=|H\rangle,$ (46)

which is equivalent to

$M_{{\cal N}^{\dagger}}^{\otimes k}|Q\rangle\\!\langle Q|M_{{\cal
N}^{\dagger}}^{\otimes k\dagger}=|H\rangle\\!\langle H|.$ (47)

Then take partial trace on both side,

$\operatorname{Tr}_{B}[M_{{\cal N}^{\dagger}}^{\otimes k}|Q\rangle\\!\langle
Q|M_{{\cal N}^{\dagger}}^{\otimes
k\dagger}]=\operatorname{Tr}_{B}[|H\rangle\\!\langle H|].$ (48)

We can simplify it into

$M_{{\cal N}^{\dagger}}\operatorname{Tr}_{B}[I\otimes M^{\otimes(k-1)}_{{\cal
N}^{\dagger}}|Q\rangle\\!\langle Q|I\otimes M_{{\cal
N}^{\dagger}}^{\otimes(k-1)\dagger}]M_{{\cal
N}^{\dagger}}^{\dagger}=\operatorname{Tr}_{B}[|H\rangle\\!\langle H|].$ (49)

For further step, let’s take $P=\operatorname{Tr}_{B}[I\otimes
M^{\otimes(k-1)}_{{\cal N}^{\dagger}}|Q\rangle\\!\langle Q|I\otimes M_{{\cal
N}^{\dagger}}^{\otimes(k-1)\dagger}]$, and
$T=\operatorname{Tr}_{B}[|H\rangle\\!\langle H|]$, thus

$M_{{\cal N}^{\dagger}}PM_{{\cal N}^{\dagger}}^{\dagger}=T.$ (50)

From the Lemma of full effective rank, the matrix $T$ is full rank, thus the
rank of $M_{{\cal N}^{\dagger}}PM_{{\cal N}^{\dagger}}^{\dagger}$ is full as
well. Due to the fact that
$\operatorname{Rank}(AB)\leq\min(\operatorname{Rank}(A),\operatorname{Rank}(B))$,
the rank on the left hand side is

$\operatorname{Rank}(M_{{\cal N}^{\dagger}}PM_{{\cal
N}^{\dagger}}^{\dagger})\leq\min\\{\operatorname{Rank}(M_{{\cal
N}^{\dagger}}),\operatorname{Rank}(P),\operatorname{Rank}(M_{{\cal
N}^{\dagger}}^{\dagger})\\},$ (51)

implying that $M_{{\cal N}^{\dagger}}$, $P$ and $M_{{\cal
N}^{\dagger}}^{\dagger}$ have full ranks. Then we have

$M_{{\cal N}^{\dagger}}=\sum_{k}\bar{E_{k}}^{\dagger}\otimes
E_{k}^{\dagger}=\sum_{k}(\bar{E_{k}}\otimes E_{k})^{\dagger}=M_{\cal
N}^{\dagger},$ (52)

since conjugate transpose won’t affect the rank of matrix, which means
$M_{\cal N}$ has full rank. This concludes that ${\cal N}$ is invertible, and
eventually an contradiction is formed. Thus, no quantum protocol exists to
extract high-order moment $\operatorname{Tr}[\rho^{k}]$.  $\sqcap$$\sqcup$

##  Appendix D Proof for Lemma 2

###### Lemma 10

Zhao  _et al._ (2023) Given a quantum channel ${\cal N}$ and an observable
$O$, there exists a HPTS map ${\cal D}$ such that ${\cal
N}^{\dagger}\circ{\cal D}^{\dagger}(O)=O$ if and only if $O\in{\cal
N}^{\dagger}({\cal L}^{H}_{A^{\prime}})$.

Lemma 2 Given an invertible quantum channel ${\cal N}$ and an observable $O$,
there exists a quantum channel ${\cal C}$ and real coefficients $t,f$ such
that

${\cal N}^{\dagger}\circ{\cal C}^{\dagger}(O)=\frac{1}{f}\left(O+tI\right).$
(53)

###### Proof.

Since the noise channel ${\cal N}$ is invertible, the image of the adjoint map
of the noise channel ${\cal N}^{\dagger}({\cal L})$ has full dimension and
hence $O\in{\cal N}^{\dagger}({\cal L})$. From Lemma 10, there exists a HPTS
map ${\cal D}$ such that ${\cal N}^{\dagger}\circ{\cal D}^{\dagger}(H)=H$.
Then we take the trick of observable shift, which is

$\displaystyle{\cal N}^{\dagger}\left({\cal D}^{\dagger}(H)\right)=H+tI-tI.$
(54)

Since ${\cal N}$ is CPTP, its adjoint map ${\cal N}^{\dagger}$ is unital-
preserving, implying

${\cal N}^{\dagger}\circ{\cal D}^{\dagger}(H+tI)=H+tI.$ (55)

Denote $\tilde{{\cal D}}^{\dagger}(H)={\cal D}^{\dagger}(H)+tI$. Since ${\cal
D}$ is an HPTS map, thus ${\cal D}^{\dagger}$ is a HPUS, i.e., ${\cal
D}^{\dagger}(I)=aI$, where $a$ is a real coefficient. Correspondingly,

$\tilde{{\cal D}}^{\dagger}(I)={\cal D}^{\dagger}(I)+tI=(a+t)I$ (56)

is also a HPUS map. Due to the fact that the adjoint of an HPUS map is HPTS,
meaning $\tilde{{\cal D}}$ is HPTS map. As long as the value $t$ no smaller
than the absolute minimum eigenvalue of ${\cal D}^{\dagger}(H)$, i.e.,
$t\geq|\min\\{eig({\cal D}^{\dagger}(H))\\}|$, the map $\tilde{{\cal D}}$
reduces to CPTS. Thus, we have

${\cal N}^{\dagger}\circ\tilde{{\cal D}}^{\dagger}(H)=H+tI$ (57)

Note that any CPTS map and be written as a coefficient $f$ times a CPTP map
${\cal C}$, which is $\tilde{{\cal D}}=f{\cal C}$, meaning

$\displaystyle f{\cal N}^{\dagger}\circ{\cal C}^{\dagger}(H)$
$\displaystyle=H+tI$ (58) $\displaystyle\Rightarrow{\cal
N}^{\dagger}\circ{\cal C}^{\dagger}(H)$ $\displaystyle=\frac{1}{f}(H+tI)$ (59)

which complete the proof.  $\sqcap$$\sqcup$

##  Appendix E Dual SDP for sampling overhead

The original SDP is

$\displaystyle f_{\min}({\cal N},k)$ $\displaystyle=\min\quad f$ (60a) subject
to $\displaystyle\quad J_{\tilde{{\cal C}}}\geq 0$ (60b)
$\displaystyle\quad\operatorname{Tr}_{C}[J_{\tilde{{\cal C}}_{BC}}]=fI_{B}$
(60c) $\displaystyle\quad J_{{\cal
F}_{AC}}\equiv\operatorname{Tr}_{B}[(J_{{\cal N}_{AB}^{\otimes
k}}^{T_{B}}\otimes I_{C})(I_{A}\otimes J_{\tilde{{\cal C}}_{BC}})]$ (60d)
$\displaystyle\quad\operatorname{Tr}_{C}[(I_{A}\otimes H_{C}^{T})J_{{\cal
F}_{AC}}^{T}]=H_{A}+tI_{A}$ (60e)

The lagrange function is

$\displaystyle L(f,J_{\tilde{{\cal C}}},t,M,K)$ $\displaystyle=f+\langle
M,\operatorname{Tr}_{C}[J_{\tilde{{\cal C}}}]-fI_{B}\rangle+\langle
K,\operatorname{Tr}_{C}[(I_{A}\otimes H_{C}^{T})J_{{\cal
F}_{AC}}^{T}]-H_{A}-tI_{A}\rangle$ (61)
$\displaystyle=f(1-\operatorname{Tr}[M])+\langle M\otimes I,J_{\tilde{{\cal
C}}}\rangle-\langle K,H\rangle-t\operatorname{Tr}[K]+\langle
K,\operatorname{Tr}_{C}[(I_{A}\otimes H_{C}^{T})J_{{\cal
F}_{AC}}^{T}]\rangle,$ (62)

where $M,K$ are the dual variables. The last term can be expanded as

$\displaystyle\langle K,\operatorname{Tr}_{A}[(H_{A}^{T}\otimes I_{C})J_{{\cal
F}_{AC}}^{T}]\rangle$
$\displaystyle=\operatorname{Tr}[K\operatorname{Tr}_{C}[(I_{A}\otimes
H_{C}^{T})J_{{\cal F}_{AC}}^{T}]]$ (63)
$\displaystyle=\operatorname{Tr}[K\otimes I(I_{A}\otimes H_{C}^{T})J_{{\cal
F}_{AC}}^{T}]$ (64) $\displaystyle=\operatorname{Tr}[(K_{A}\otimes
H_{C}^{T})J_{{\cal F}_{AC}}^{T}]$ (65)
$\displaystyle=\operatorname{Tr}[(K_{A}^{T}\otimes H_{C})J_{{\cal F}_{AC}}]$
(66) $\displaystyle=\operatorname{Tr}[(K_{A}^{T}\otimes
H_{C})\operatorname{Tr}_{B}[(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes
I_{C})(I_{A}\otimes J_{\tilde{{\cal C}}_{BC}})]]$ (67)
$\displaystyle=\operatorname{Tr}[(K_{A}^{T}\otimes I_{B}\otimes
H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})(I_{A}\otimes
J_{\tilde{{\cal C}}_{BC}})]$ (68)
$\displaystyle=\operatorname{Tr}[\operatorname{Tr}_{A}[(K_{A}^{T}\otimes
I_{B}\otimes H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes
I_{C})]J_{\tilde{{\cal C}}_{BC}}]$ (69)

Thus, the Lagrange function can be expressed as

$\displaystyle{\cal L}(f,J_{\tilde{{\cal C}}},t,M,K)$
$\displaystyle=f(1-\operatorname{Tr}[M])+\langle M\otimes I,J_{\tilde{{\cal
C}}}\rangle-\langle
K,H\rangle-t\operatorname{Tr}[K]+\langle\operatorname{Tr}_{A}[(K_{A}^{T}\otimes
I_{B}\otimes H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes
I_{C})],J_{\tilde{{\cal C}}_{BC}}\rangle$ (70)
$\displaystyle=f(1-\operatorname{Tr}[M])-\langle
K,H\rangle-t\operatorname{Tr}[K]+\langle M\otimes
I+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes H_{C})(J_{{\cal
N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})],J_{\tilde{{\cal C}}}\rangle$ (71)

The corresponding Lagrange dual function is

$\displaystyle f(M,K)$ $\displaystyle=\inf_{J_{\tilde{{\cal D}}},t,f}{\cal
L}(f,J_{\tilde{{\cal C}}},t,M,K)$ (72) $\displaystyle=-\langle
K,H\rangle+\inf_{J_{\tilde{{\cal
D}}},t,f}(f(1-\operatorname{Tr}[M])-t\operatorname{Tr}[K]+\langle M\otimes
I+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes H_{C})(J_{{\cal
N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})],J_{\tilde{{\cal C}}}\rangle)$ (73)

In order to bound the infimum term, every sub-term inside the infimum should
greater greater or equal to 0, thus, we have

$\displaystyle 1-\operatorname{Tr}[M]\geq 0,\quad t\operatorname{Tr}[K]=0$
(74) $\displaystyle\langle M\otimes I+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes
I_{B}\otimes H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes
I_{C})],J_{\tilde{{\cal C}}}\rangle\geq 0$ (75)

Since $J_{\tilde{{\cal C}}}\geq 0$, we have $M\otimes
I+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes H_{C})(J_{{\cal
N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})]\geq 0$. Thus, we arrive at the
following dual SDP:

$\displaystyle f_{\min}({\cal N},k)$
$\displaystyle=\max\quad-\operatorname{Tr}[KH]$ (76a) subject to
$\displaystyle\quad\operatorname{Tr}[M]\leq 1$ (76b)
$\displaystyle\quad\operatorname{Tr}[K]=0$ (76c) $\displaystyle\quad M\otimes
I+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes H_{C})(J_{{\cal
N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})]\geq 0$ (76d)

##  Appendix F Optimal Sampling Cost for depolarizing channel (Proposition 4)

The depolarizing channel is one the most common noise in quantum computers
which maps part of quantum state into maximally mixed state. The $1$-qubit
depolarizing channel can be written in the following form:

${\cal N}^{\epsilon}_{\rm DE}(\rho)=(1-\epsilon)\rho+\epsilon\frac{I}{2},$
(77)

where $\epsilon\in[0,1]$ is the noisy level and $I$ is identity matrix. The
retrieving protocol and retrieving cost for such noise channel was well
studied in Zhao  _et al._ (2023). Here we will focus on retrieving the second
order information $\operatorname{Tr}[\rho^{2}]$ from two-copy of depolarizing
noised state ${\cal N}(\rho\otimes\rho)={\cal N}^{\epsilon}_{\rm
DE}\otimes{\cal N}^{\epsilon}_{\rm DE}(\rho\otimes\rho)$, after applying such
channel, the state becomes

$\displaystyle{\cal N}(\rho\otimes\rho)$ $\displaystyle={\cal
N}^{\epsilon}_{\rm DE}(\rho)\otimes{\cal N}^{\epsilon}_{\rm DE}(\rho)$ (78)
$\displaystyle=[(1-\epsilon)\rho+\epsilon\frac{I}{2}]\otimes[(1-\epsilon)\rho+\epsilon\frac{I}{2}]$
(79)
$\displaystyle=(1-\epsilon)^{2}\rho\otimes\rho+\frac{\epsilon(1-\epsilon)}{2}I\otimes\rho+\frac{(1-\epsilon)\epsilon}{2}\rho\otimes
I+\frac{\epsilon^{2}}{4}I\otimes I.$ (80)

Proposition 4 Given noisy states ${\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes
2}$, and error tolerance $\delta$, the second order moment
$\operatorname{Tr}[\rho^{2}]$ can be estimated by $f\operatorname{Tr}[H{\cal
C}\circ{\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes 2}]-t$, with optimal sample
complexity $\mathcal{O}(1/(\delta^{2}(1-\epsilon)^{4})$, where
$f=\frac{1}{(1-\epsilon)^{2}}$,
$t=\frac{1-(1-\epsilon)^{2}}{2(1-\epsilon)^{2}}$. The term
$\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes
2}]$ can be estimated by implementing a quantum retriever ${\cal C}$ on noisy
states and making measurement over moment observable $H$. Moreover, there
exists an ensemble of unitary operations $\\{p_{j},U_{j}\\}_{j}$ such that the
action of the retriever ${\cal C}$ can be interpreted as,

${\cal C}(\cdot)=\sum_{j=1}p_{j}U_{j}\cdot U_{j}^{\dagger}$ (81)

###### Proof.

To begin with, we give the explicit form of the retriever ${\cal C}$.

${\cal C}(\cdot)=\sum_{j=1}^{12}\frac{1}{12}U_{j}\cdot U_{j}^{\dagger},$ (82)

where all the probabilities $p_{j}$ are the same, i.e., $p_{j}=\frac{1}{12}$,
and the corresponding unitaries $U_{j}$ are

$\displaystyle U_{1}$ $\displaystyle=I\otimes I;\quad U_{2}=X\otimes X;\quad
U_{3}=Y\otimes Y;\quad U_{4}=Z\otimes Z;$ (83) $\displaystyle U_{5}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}-i&1&1&i\\\ i&1&-1&i\\\ i&-1&1&i\\\
-i&-1&-1&i\end{pmatrix};\ U_{6}=\frac{1}{2}\begin{pmatrix}i&-1&-1&-i\\\
i&1&-1&i\\\ i&-1&1&i\\\ i&1&1&-i\end{pmatrix};\
U_{7}=\frac{1}{2}\begin{pmatrix}i&i&i&i\\\ -1&1&-1&1\\\ -1&-1&1&1\\\
-i&i&i&-i\end{pmatrix};U_{8}=\frac{1}{2}\begin{pmatrix}-i&i&i&-i\\\
1&1&-1&-1\\\ 1&-1&1&-1\\\ i&i&i&i\end{pmatrix};$ $\displaystyle U_{9}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}i&1&1&-i\\\ -i&1&-1&-i\\\
-i&-1&1&-i\\\ i&-1&-1&-i\end{pmatrix};\
U_{10}=\frac{1}{2}\begin{pmatrix}-i&-1&-1&i\\\ -i&1&-1&-i\\\ -i&-1&1&-i\\\
-i&1&1&i\end{pmatrix};U_{11}=\frac{1}{2}\begin{pmatrix}i&-i&-i&i\\\
1&1&-1&-1\\\ 1&-1&1&-1\\\ -i&-i&-i&-i\end{pmatrix};\
U_{12}=\frac{1}{2}\begin{pmatrix}-i&-i&-i&-i\\\ -1&1&-1&1\\\ -1&-1&1&1\\\
i&-i&-i&i\end{pmatrix}.$

Next, we are going to prove the retriever ${\cal C}$ has the optimal sampling
overhead. The optimal complexity $\mathcal{O}(1/(\delta^{2}(1-\epsilon)^{4})$
implies the optimal sampling overhead is $\frac{1}{(1-\epsilon)^{2}}$. To
prove estimating the second order moment by implementing quantum channel $C$
with optimal sampling overhead is $\frac{1}{(1-\epsilon)^{2}}$, we divide the
process in to two step. For the first step, we prove that quantum channel
${\cal C}$ is a feasible solution to the prime problem with sampling overhead
$\frac{1}{(1-\epsilon)^{2}}$, meaning the optimal overhead $\gamma^{*}$,
$\gamma^{*}\leq\frac{1}{(1-\epsilon)^{2}}$. Then we are going to find one
feasible solution to dual problem with overhead $\frac{1}{(1-\epsilon)^{2}}$,
meaning $\gamma^{*}\geq\frac{1}{(1-\epsilon)^{2}}$. Thus we have the optimal
overhead $\gamma^{*}=\frac{1}{(1-\epsilon)^{2}}$, which complete the proof.

Now, we are going proof the first part. At first, we need to write the Choi
matrix form of the retriever ${\cal C}$, which is

$J_{\cal C}=\frac{1}{4}IIII+\frac{1}{12}(XX+YY+ZZ)\otimes(XX+YY+ZZ).$ (84)

For an arbitrary state $\rho\otimes\rho$, after applying the noise ${\cal N}$
onto it, we denote the state as $\rho^{\prime}={\cal N}(\rho\otimes\rho)$,
thus

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$ $\displaystyle={\cal
C}(\rho^{\prime})=\operatorname{Tr}_{A}[(\rho^{\prime T}\otimes II)J_{\cal
C}]$ (85) $\displaystyle=\operatorname{Tr}_{A}\Big{[}(\rho^{\prime T}\otimes
II)(\frac{1}{4}IIII+\frac{1}{12}(XX+YY+ZZ)\otimes(XX+YY+ZZ))$ (86)
$\displaystyle=\frac{1}{4}\operatorname{Tr}[\rho^{\prime
T}]II+\frac{1}{12}\operatorname{Tr}[\rho^{\prime T}(XX+YY+ZZ)](XX+YY+ZZ))$
(87) $\displaystyle=\frac{1}{4}II+\frac{1}{12}\operatorname{Tr}[\rho^{\prime
T}(XX+YY+ZZ)](XX+YY+ZZ))\Big{]}$ (88)

Note that above equations utilized the face that transpose operation is trace
preserving, i.e., $\operatorname{Tr}[\rho^{T}]=\operatorname{Tr}[\rho]=1$.
Since the matrix $XX+YY+ZZ$ is symmetry, thus we have

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$
$\displaystyle=\frac{1}{4}II+\frac{1}{12}\operatorname{Tr}[\rho^{\prime
T}(XX+YY+ZZ)^{T}](XX+YY+ZZ))$ (89)
$\displaystyle=\frac{1}{4}II+\frac{1}{12}\operatorname{Tr}[\rho^{\prime}(XX+YY+ZZ)](XX+YY+ZZ))$
(90)

The trace term can be calculated by substituting Eq. (80), we have

$\displaystyle\operatorname{Tr}[\rho^{\prime}(XX+YY+ZZ)]$
$\displaystyle=\operatorname{Tr}[((1-\epsilon)^{2}\rho\otimes\rho+\frac{\epsilon(1-\epsilon)}{2}I\otimes\rho+\frac{(1-\epsilon)\epsilon}{2}\rho\otimes
I+\epsilon^{2}(I\otimes I))(XX+YY+ZZ)]$ (91)
$\displaystyle=(1-\epsilon)^{2}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)]+\frac{\epsilon(1-\epsilon)}{2}\operatorname{Tr}[X\otimes\rho
X+Y\otimes\rho Y+Z\otimes\rho Z]$
$\displaystyle\quad+\frac{(1-\epsilon)\epsilon}{2}\operatorname{Tr}[\rho
X\otimes X+\rho Y\otimes Y+\rho Z\otimes
Z]+\frac{\epsilon^{2}}{4}\operatorname{Tr}[XX+YY+ZZ]$ (92)
$\displaystyle=(1-\epsilon)^{2}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)].$
(93)

In the second equation, since $X,Y,Z$ are all traceless Hermitian matrices,
all terms are zeros except the first term. Replace the equation back to Eq.
(90), then

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$
$\displaystyle=\frac{1}{4}II+\frac{(1-\epsilon)^{2}}{12}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)](XX+YY+ZZ).$
(94)

The information $\operatorname{Tr}[\rho^{2}]$ is estimated from
$\operatorname{Tr}[H(\rho\otimes\rho)]$, where $H=\frac{1}{2}(II+XX+YY+ZZ)$ is
cyclic permutation operator (in the 2-qubit case, $H$ is just SWAP operator).
It is easy to check that

$\displaystyle\operatorname{Tr}[H{\cal C}\circ{\cal N}(\rho\otimes\rho)]$
$\displaystyle=\frac{1}{(1-\epsilon)^{2}}\Big{[}\frac{1}{4}\operatorname{Tr}[H*II]+\frac{(1-\epsilon)^{2}}{12}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)]\operatorname{Tr}[H*(XX+YY+ZZ)]\Big{]}.$
(95)

We can quickly get $\operatorname{Tr}[H]=2$ and
$\operatorname{Tr}[H*(XX+YY+ZZ)]=6$. Then

$\displaystyle\operatorname{Tr}[H{\cal C}\circ{\cal N}(\rho\otimes\rho)]$
$\displaystyle=\frac{1}{2}+\frac{(1-\epsilon)^{2}}{2}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)]$
(96)
$\displaystyle=(1-\epsilon)^{2}\Big{[}\frac{1}{2(1-\epsilon)^{2}}+\frac{1}{2}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)]\Big{]}$
(97)
$\displaystyle=(1-\epsilon)^{2}\Big{[}\frac{1}{2(1-\epsilon)^{2}}-\frac{1}{2}+\frac{1}{2}+\frac{1}{2}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ)]\Big{]}$
(98)
$\displaystyle=(1-\epsilon)^{2}\Big{[}\frac{1}{2(1-\epsilon)^{2}}-\frac{(1-\epsilon)^{2}}{2(1-\epsilon)^{2}}+\frac{1}{2}\operatorname{Tr}[\rho\otimes\rho(XX+YY+ZZ+II)]\Big{]}$
(99)
$\displaystyle=(1-\epsilon)^{2}\Big{[}\frac{2\epsilon^{2}-\epsilon^{2}}{2(1-\epsilon)^{2}}+\operatorname{Tr}[\rho\otimes\rho
H]\Big{]}$ (100)
$\displaystyle=(1-\epsilon)^{2}\Big{[}\frac{2\epsilon^{2}-\epsilon^{2}}{2(1-\epsilon)^{2}}+\operatorname{Tr}[\rho^{2}]\Big{]}.$
(101)

The desired high-order moment $\operatorname{Tr}[\rho^{2}]$ value equals to

$\operatorname{Tr}[\rho^{2}]=f\operatorname{Tr}[H{\cal C}\circ{\cal
N}(\rho\otimes\rho)]-t$ (102)

where $f=\frac{1}{(1-\epsilon)^{2}}$ is the sampling overhead and
$t=\frac{1-(1-\epsilon)^{2}}{2(1-\epsilon)^{2}}$ is the shifted distance. In
order to estimate the value with the error $\delta$, the sampling overhead
should be $1/(1-\epsilon)^{2}$. Thus, we have $\gamma^{*}\leq
1/(1-\epsilon)^{2}$.

Next, we are going to use dual SDP to show that $\gamma^{*}\geq
1/(1-\epsilon)^{2}$. We set the dual variables as
$M=\frac{1}{4}II-\frac{1}{12}(XX+YY+ZZ)$, and $K=q(XX+YY+ZZ)$ where
$q=-\frac{1}{6(1-\epsilon)^{2}}$. We will show the variables $\\{M,K\\}$ is a
feasible solution to the dual problem.

If substitute the variables into the dual problem Eq. (76), we can easily
check that $\operatorname{Tr}[M]\leq 1,\operatorname{Tr}[K]=0$. For the last
condition, we have

$\displaystyle\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes
H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes
I_{C})]=-\frac{1}{12}(XX+YY+ZZ)(II+XX+YY+ZZ)$ (103)

Thus, we have

$\displaystyle M\otimes II+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes
I_{B}\otimes H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})]$
(104)
$\displaystyle=\frac{1}{4}IIII-\frac{1}{6}(XX+YY+ZZ)II-\frac{1}{12}(XX+YY+ZZ)(XX+YY+ZZ)$
(105) $\displaystyle\geq 0,$ (106)

which means $\\{M,K\\}$ is a feasible solution to the dual SDP. Therefore, we
have $\gamma*\geq-\operatorname{Tr}[KH]=\frac{1}{(1-\epsilon)^{2}}$. Combined
with prime part, we have $\gamma^{*}=\frac{1}{(1-\epsilon)^{2}}$, meaning that
given error tolerance $\delta$, the optimal sample compleixty is ${\cal
O}(\frac{1}{\delta^{2}(1-\epsilon)^{4}})$. The proof is complete.
$\sqcap$$\sqcup$

##  Appendix G Optimal Sampling cost for amplitude damping channel
(Proposition 5)

The amplitude damping (AD) channel is a physical channel that describes the
energy leakage, dropping from a high energy state to a low energy state. The
qubit amplitude damping channel ${\cal N}_{\rm AD}^{\varepsilon}$ has two
Kraus operators: $A_{0}^{\varepsilon}\coloneqq|0\rangle\\!\langle
0|+\sqrt{1-\varepsilon}|1\rangle\\!\langle 1|$ and
$A_{1}^{\varepsilon}\coloneqq\sqrt{\varepsilon}|0\rangle\\!\langle 1|$. A
single qubit state $\begin{pmatrix}\rho_{00}&\rho_{01}\\\
\rho_{10}&\rho_{11}\end{pmatrix}$ after going through the AD channel is

${\cal N}_{\rm
AD}^{\varepsilon}(\rho)=\begin{pmatrix}\rho_{00}+\varepsilon\rho_{11}&\sqrt{1-\varepsilon}\rho_{01}\\\
\sqrt{1-\varepsilon}\rho_{10}&(1-\varepsilon)\rho_{11}\end{pmatrix}$ (107)

where $\varepsilon\in[0,1]$ is the damping factor. This part we will focus on
retrieving the information $\operatorname{Tr}[\rho^{2}]$ from two-copy AD
channel noised state ${\cal N}_{\rm AD}^{\varepsilon}\otimes{\cal N}_{\rm
AD}^{\varepsilon}(\rho\otimes\rho)$. The noised state is

$\displaystyle{\cal N}(\rho\otimes\rho)$ $\displaystyle={\cal N}_{\rm
AD}^{\varepsilon}\otimes{\cal N}_{\rm AD}^{\varepsilon}(\rho\otimes\rho)$
(108)
$\displaystyle=\begin{pmatrix}(\rho_{00}+\varepsilon\rho_{11})^{2}&(\rho_{00}+\varepsilon\rho_{11})(\sqrt{1-\varepsilon}\rho_{01})&(\sqrt{1-\varepsilon}\rho_{01})(\rho_{00}+\varepsilon\rho_{11})&(1-\varepsilon)\rho_{01}^{2}\\\
(\rho_{00}+\varepsilon\rho_{11})(\sqrt{1-\varepsilon}\rho_{10})&(\rho_{00}+\varepsilon\rho_{11})(1-\varepsilon)\rho_{11}&(1-\varepsilon)\rho_{01}\rho_{10}&(1-\varepsilon)\sqrt{1-\varepsilon}\rho_{01}\rho_{11}\\\
\sqrt{1-\varepsilon}\rho_{10}(\rho_{00}+\varepsilon\rho_{11})&(1-\varepsilon)\rho_{10}\rho_{01}&(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})&(1-\varepsilon)\sqrt{1-\varepsilon}\rho_{11}\rho_{01}\\\
(1-\varepsilon)\rho_{10}^{2}&(1-\varepsilon)\sqrt{1-\varepsilon}\rho_{10}\rho_{11}&(1-\varepsilon)\sqrt{1-\varepsilon}\rho_{11}\rho_{10}&(1-\varepsilon)^{2}\rho_{11}^{2}\\\
\end{pmatrix}$ (109)

Proposition 5 Given noisy states ${\cal N}^{\varepsilon}_{\rm
AD}(\rho)^{\otimes 2}$, and error tolerance $\delta$, the second order moment
$\operatorname{Tr}[\rho^{2}]$ can be estimated by $f\operatorname{Tr}[H{\cal
C}\circ{\cal N}^{\varepsilon}_{\rm AD}(\rho)^{\otimes 2}]-t$, with optimal
sample complexity $\mathcal{O}(1/(\delta^{2}(1-\varepsilon)^{4})$, where
$f=\frac{1}{(1-\varepsilon)^{2}}$,
$t=-\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}$. The term
$\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\varepsilon}_{\rm
AD}(\rho)^{\otimes 2}]$ can be estimated by implementing a quantum retriever
${\cal C}$ on noisy states and making measurement. Moreover, the Choi matrix
of such the retriever ${\cal C}$ is

$\displaystyle J_{\cal C}$ $\displaystyle=|00\rangle\\!\langle
00|\otimes\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$
$\displaystyle\quad+|\Psi^{+}\rangle\\!\langle\Psi^{+}|\otimes\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$
$\displaystyle\quad+|\Psi^{-}\rangle\\!\langle\Psi^{-}|\otimes\frac{1}{2}(II-H)$
$\displaystyle\quad+|11\rangle\\!\langle 11|\otimes\frac{1}{6}(II+H).$ (110)

where
$|\Psi^{\pm}\rangle\\!\langle\Psi^{\pm}|=\frac{1}{2}(|01\rangle\pm|10\rangle)(\langle
01|\pm\langle 10|)$ are Bell states.

###### Proof.

The optimal complexity $\mathcal{O}(1/(\delta^{2}(1-\varepsilon)^{4})$ implies
the optimal sampling overhead is $\frac{1}{(1-\varepsilon)^{2}}$. To prove to
estimate the second order moment by implementing quantum channel $C$ with
optimal sampling overhead is $\frac{1}{(1-\varepsilon)^{2}}$, we divide the
process into two steps. For the first step, we prove that quantum channel
${\cal C}$ is a feasible solution to the prime problem with sampling overhead
$\frac{1}{(1-\varepsilon)^{2}}$, meaning the optimal overhead $\gamma^{*}$,
$\gamma^{*}\leq\frac{1}{(1-\varepsilon)^{2}}$. Then we are going to find one
feasible solution to the dual problem with overhead
$\frac{1}{(1-\varepsilon)^{2}}$, meaning
$\gamma^{*}\geq\frac{1}{(1-\varepsilon)^{2}}$. Thus we have the optimal
overhead $\gamma^{*}=\frac{1}{(1-\varepsilon)^{2}}$, which complete the proof.

Now, we are going to prove the first part. For an arbitrary state
$\rho\otimes\rho$, after applying the noise ${\cal N}$ onto it, we denote the
state as $\rho^{\prime}={\cal N}(\rho\otimes\rho)$, thus

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$ $\displaystyle={\cal
C}(\rho^{\prime})=\operatorname{Tr}_{A}[(\rho^{\prime T}\otimes II)J_{\cal
C}]$ $\displaystyle=\operatorname{Tr}_{A}\Big{[}(\rho^{\prime T}\otimes
II)\Big{(}|00\rangle\\!\langle
00|\otimes\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$
$\displaystyle\quad+\frac{1}{2}(|01\rangle+|10\rangle)(\langle 01|+\langle
10|)\otimes\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)$
$\displaystyle\quad+\frac{1}{2}(|01\rangle-|10\rangle)(\langle 01|-\langle
10|)\otimes\frac{1}{2}(II-H)$ $\displaystyle\quad+|11\rangle\\!\langle
11|\otimes\frac{1}{6}(II+H)\Big{)}\Big{]}$ (111)
$\displaystyle=\operatorname{Tr}[\rho^{\prime T}|00\rangle\\!\langle
00|]\bigg{(}\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)\bigg{)}$
$\displaystyle\quad+\operatorname{Tr}[\rho^{\prime
T}((|01\rangle+|10\rangle)(\langle 01|+\langle
10|)]\bigg{(}\frac{1}{12}((1+2\varepsilon)II+(1-4\varepsilon)H)\bigg{)}$
$\displaystyle\quad+\operatorname{Tr}[\rho^{\prime
T}((|01\rangle-|10\rangle)(\langle 01|-\langle
10|)]\bigg{(}\frac{1}{4}(II-H)\bigg{)}$
$\displaystyle\quad+\operatorname{Tr}[\rho^{\prime T}|11\rangle\\!\langle
11|]\bigg{(}\frac{1}{6}(II+H)\bigg{)}$ (112)

we can take transpose of the quantum noisy state $\rho^{\prime}$, as shown in
Eq. (109), then substitute it into our equation, we get

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$
$\displaystyle=(\rho_{00}+\varepsilon\rho_{11})^{2}\bigg{(}\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)\bigg{)}$
$\displaystyle\quad+[(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})+(1-\varepsilon)\rho_{01}\rho_{10}]\bigg{(}\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H)\bigg{)}$
$\displaystyle\quad+[(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})-(1-\varepsilon)\rho_{01}\rho_{10}]\bigg{(}\frac{1}{2}(II-H)\bigg{)}$
$\displaystyle\quad+(1-\varepsilon)^{2}\rho_{11}^{2}\bigg{(}\frac{1}{6}(II+H)\bigg{)}$
(113)

The value $\operatorname{Tr}[\rho^{2}]$ is usually estimated by
$\operatorname{Tr}[H(\rho\otimes\rho)]$, where $H=\frac{1}{2}(II+XX+YY+ZZ)$ is
the two-qubit cyclic permutation operator. Then the estimated value from our
protocol can be arrived at

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$
$\displaystyle=(\rho_{00}+\varepsilon\rho_{11})^{2}\bigg{(}\frac{1}{6}\operatorname{Tr}[H((1+2\varepsilon)II+(1-4\varepsilon)H)]\bigg{)}$
$\displaystyle\quad+[(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})+(1-\varepsilon)\rho_{01}\rho_{10}]\bigg{(}\frac{1}{6}\operatorname{Tr}[H((1+2\varepsilon)II+(1-4\varepsilon)H)]\bigg{)}$
$\displaystyle\quad+[(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})-(1-\varepsilon)\rho_{01}\rho_{10}]\bigg{(}\frac{1}{2}\operatorname{Tr}[H(II-H)]\bigg{)}$
$\displaystyle\quad+(1-\varepsilon)^{2}\rho_{11}^{2}\bigg{(}\frac{1}{6}\operatorname{Tr}[H(II+H)]\bigg{)}$
(114)

Note that the $\operatorname{Tr}[H]=2$ and $\operatorname{Tr}[H\cdot H]=4$,
then we have

$\displaystyle\operatorname{Tr}[H{\cal C}\circ{\cal N}(\rho\otimes\rho)]$
$\displaystyle=(\rho_{00}+\varepsilon\rho_{11})^{2}(1-2\varepsilon))+[(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})+(1-\varepsilon)\rho_{01}\rho_{10}](1-2\varepsilon)$
$\displaystyle\quad-[(1-\varepsilon)\rho_{11}(\rho_{00}+\varepsilon\rho_{11})-(1-\varepsilon)\rho_{01}\rho_{10}]+(1-\varepsilon)^{2}\rho_{11}^{2}$
(115)
$\displaystyle=(1-2\varepsilon)\rho_{00}^{2}+2(1-\varepsilon)^{2}\rho_{01}\rho_{10}-2\varepsilon^{2}\rho_{00}\rho_{11}+(1-2\varepsilon)\rho_{11}^{2}$
(116)
$\displaystyle=(1-\varepsilon)^{2}\rho_{00}^{2}+2(1-\varepsilon)^{2}\rho_{01}\rho_{10}+(1-\varepsilon)^{2}\rho_{11}^{2}-\varepsilon^{2}\rho_{00}^{2}+2\varepsilon^{2}\rho_{00}\rho_{11}+\varepsilon^{2}\rho_{11}^{2}$
(117)
$\displaystyle=(1-\varepsilon)^{2}\Big{[}\rho_{00}^{2}+\rho_{11}^{2}+2\rho_{10}\rho_{01}-\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}[\rho_{00}^{2}+\rho_{11}^{2}+\rho_{00}\rho_{11}+\rho_{11}\rho_{00}]\Big{]}$
(118)
$\displaystyle=(1-\varepsilon)^{2}\Big{[}\operatorname{Tr}[H\rho\otimes\rho]-\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}\Big{]}$
(119)
$\displaystyle=(1-\varepsilon)^{2}\Big{[}\operatorname{Tr}[\rho^{2}]-\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}\Big{]}.$
(120)

The desired high-order moment $\operatorname{Tr}[\rho^{2}]$ value equals to

$\operatorname{Tr}[\rho^{2}]=f\operatorname{Tr}[H{\cal C}\circ{\cal
N}(\rho\otimes\rho)]-t$ (122)

where $f=\frac{1}{(1-\varepsilon)^{2}}$ is the sampling overhead and
$t=-\frac{\varepsilon^{2}}{(1-\varepsilon)^{2}}$ is the shifted distance. In
order to estimate the value with the error $\delta$, the sampling overhead
should be $1/(1-\varepsilon)^{2}$. Thus, we have $\gamma^{*}\leq
1/(1-\varepsilon)^{2}$

Next, we are going to use dual SDP to show that $\gamma^{*}\geq
1/(1-\varepsilon)^{2}$. We set the dual variables as

$\displaystyle M$ $\displaystyle=\frac{1}{4}(|01\rangle-|10\rangle)(\langle
01|-\langle 10|)+\frac{1}{2}|11\rangle\\!\langle 11|$ (123) $\displaystyle K$
$\displaystyle=\frac{1}{2(1-\varepsilon)^{2}}\Big{[}-\varepsilon|00\rangle\\!\langle
00|-|11\rangle\\!\langle 11|+\frac{1+\varepsilon}{2}(|01\rangle\\!\langle
01|+|10\rangle\\!\langle 10|)+\frac{\varepsilon-1}{2}(|01\rangle\\!\langle
10|+|10\rangle\\!\langle 01|)\Big{]}.$ (124)

We will show the variables $\\{M,K\\}$ is a feasible solution to dual problem.
If substitute the variables into the dual problem Eq. (76), we can easily
check that $\operatorname{Tr}[M]\leq 1$ and $\operatorname{Tr}[K]=1$. For the
last condition, after simplifying, we have

$\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes H_{C})(J_{{\cal
N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})]=(M-|11\rangle\\!\langle
11|)\otimes H$ (125)

Thus

$\displaystyle M\otimes I+\operatorname{Tr}_{A}[(K_{A}^{T}\otimes I_{B}\otimes
H_{C})(J_{{\cal N}_{AB}^{\otimes k}}^{T_{B}}\otimes I_{C})]$ (126)
$\displaystyle=$ $\displaystyle M\otimes I+(M-|11\rangle\\!\langle 11|)\otimes
H$ (127) $\displaystyle\geq$ $\displaystyle 0$ (128)

which means $\\{M,K\\}$ is a feasible solution to the dual SDP, therefore, we
have $\gamma^{*}\geq-\operatorname{Tr}[KH]=\frac{1}{(1-\varepsilon)^{2}}$.
Combined with prime part, we have $\gamma^{*}=\frac{1}{(1-\varepsilon)^{2}}$,
meaning that given error tolerance $\delta$, the optimal sample compleixty is
${\cal O}(\frac{1}{\delta^{2}(1-\varepsilon)^{2}})$. The proof is complete.
$\sqcap$$\sqcup$

Notice in Eq. (G), for any input state $\rho^{\prime}$ to the AD-decoder
${\cal C}$, the action can be implemented as a construction of a state,

$\displaystyle{\cal C}(\rho^{\prime})$
$\displaystyle=\frac{p_{1}}{6}\left((1+2\varepsilon)II+(1-4\varepsilon)H\right)+\frac{p_{2}}{6}\left((1+2\varepsilon)II+(1-4\varepsilon)H\right)+\frac{p_{3}}{2}(II-H)+\frac{p_{4}}{6}(II+H),$
(129)

where $p_{1,2,3,4}$ can be considered as the measurement probabilities that
state $\rho^{\prime}$ could project onto the orthonormal basis states
$|00\rangle,|11\rangle,|\Psi^{+}\rangle,|\Psi^{-}\rangle$, respectively, where
$|\Psi^{\pm}\rangle=1/\sqrt{2}(|01\rangle\pm|10\rangle)$, and the four states
form an orthonormal basis ${\cal B}$. Since ${\cal C}$ is CPTP for
$\varepsilon\in[0,1]$. The output state can then be derived as a convex
combination of states,

$\sigma_{1}=\sigma_{2}=\frac{1}{6}((1+2\varepsilon)II+(1-4\varepsilon)H);\
\sigma_{3}=\frac{1}{2}(II-H);\ \sigma_{4}=\frac{1}{6}(II+H),$ (130)

which is just,

${\cal C}(\rho^{\prime})=\sum_{j=1}^{4}p_{j}\sigma_{j}.$ (131)

Now the question becomes how to construct this ensemble of states
$\sigma_{j}$’s based on ${\cal C}(\rho^{\prime})$. Notice that for states
$\sigma_{1,2}$, we have,

$\sigma_{1}=\sigma_{2}=\frac{1-\varepsilon}{3}|00\rangle\\!\langle
00|+\frac{1-\varepsilon}{3}|11\rangle\\!\langle
11|+\frac{1-\varepsilon}{3}|\Psi^{+}\rangle\\!\langle\Psi^{+}|+\varepsilon|\Psi^{-}\rangle\\!\langle\Psi^{-}|,$
(132)

and,

$\sigma_{3}=|\Psi^{-}\rangle\\!\langle\Psi^{-}|;\
\sigma_{4}=\frac{1}{3}|00\rangle\\!\langle 00|+\frac{1}{3}|11\rangle\\!\langle
11|+\frac{1}{3}|\Psi^{+}\rangle\\!\langle\Psi^{+}|.$ (133)

Notice that for these four states, the corresponding expectation value with
respect to the SWAP observable can be directly computed as

$\operatorname{Tr}[\sigma_{1,2}H]=1-2\varepsilon;\
\operatorname{Tr}[\sigma_{3}H]=-1;\ \operatorname{Tr}[\sigma_{4}H]=1.$ (134)

Figure 6: The circuit protocol for retrieving the 2nd moment
$\operatorname{Tr}[\rho^{2}]$ from AD product channel with identical noise
factor $\varepsilon$. The input state $\rho^{\prime}={\cal
N}_{\varepsilon}(\rho\otimes\rho)$ is the noise resource which is measured in
the basis ${\cal
B}=\\{|00\rangle,|11\rangle,|\Psi^{+}\rangle,|\Psi^{-}\rangle\\}$. Here we
need to include the ancillary system in the bottom 2-qubit system and use
measurement-controlled operations shown in the diagram. The middle circuit of
2-qubit then can be used to retrieve the purity $\operatorname{Tr}[\rho^{2}]$.

Notice that the four measurements, in fact, form orthogonal projective
measurements. The retrieving of $\operatorname{Tr}[\rho^{2}]$ can then be
performed via the probabilistic weighted sum of these four values based on the
measurement outcomes or equivalently realized by a circuit framework as shown
in Fig. 6. In this example, the protocol can be implemented by setting up an,
in total, $6$-qubit system where the top two are the input subsystems, the
middle two as the retrieved output subsystems and the bottom two as the
ancillary systems. We first require a $2$-qubit unitary $U$ mapping the
initial states $|00\rangle$ in the middle to an intermediate state,

$U|00\rangle=\sqrt{\frac{1}{3}}|00\rangle+\sqrt{\frac{1}{3}}|11\rangle+\sqrt{\frac{1}{3}}|\Psi^{+}\rangle.$
(135)

Later, a parameterised unitary $U_{\varepsilon}$ based on the noise factor
$\varepsilon$ is settled to create the correct state coefficients, i.e.,

$U_{\varepsilon}U|00\rangle=\sqrt{\frac{1-\varepsilon}{3}}|00\rangle+\sqrt{\frac{1-\varepsilon}{3}}|11\rangle+\sqrt{\frac{1-\varepsilon}{3}}|\Psi^{+}\rangle+\sqrt{\varepsilon}|\Psi^{-}\rangle.$
(136)

By writing out $U_{\varepsilon}$ with respect to the basis ${\cal B}$, we
could derive the following exact form as,

$\displaystyle U_{\varepsilon}=\begin{pmatrix}1&0&0&0\\\
0&\frac{3+\sqrt{3}}{6}&\frac{3-\sqrt{3}}{6}&\frac{\sqrt{3}}{3}\\\
0&\frac{3-\sqrt{3}}{6}&\frac{3+\sqrt{3}}{6}&-\frac{\sqrt{3}}{3}\\\
0&-\frac{\sqrt{3}}{3}&\frac{\sqrt{3}}{3}&\frac{\sqrt{3}}{3}\end{pmatrix}\begin{pmatrix}\sqrt{1-\varepsilon}&0&0&-\sqrt{\varepsilon}\\\
0&\sqrt{1-\varepsilon}&-\sqrt{\varepsilon}&0\\\
0&\sqrt{\varepsilon}&\sqrt{1-\varepsilon}&0\\\
\sqrt{\varepsilon}&0&0&\sqrt{1-\varepsilon}\end{pmatrix},$ (137)

where in this ${\cal B}$, the basis states are represented as,

$|00\rangle=\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\end{pmatrix};\
|11\rangle=\begin{pmatrix}0\\\ 1\\\ 0\\\ 0\end{pmatrix};\
|\Psi^{+}\rangle=\begin{pmatrix}0\\\ 0\\\ 1\\\ 0\end{pmatrix};\
|\Psi^{-}\rangle=\begin{pmatrix}0\\\ 0\\\ 0\\\ 1\end{pmatrix}.$ (138)

After matching the coefficients, a specialized control-$V$ ($CV$) gate will be
applied to generate the purification of (136) as,

$|\psi^{\prime}\rangle=CV(U_{\varepsilon}U\otimes
I)|00\rangle|00\rangle=\sqrt{\frac{1-\varepsilon}{3}}|00\rangle|00\rangle+\sqrt{\frac{1-\varepsilon}{3}}|11\rangle|01\rangle+\sqrt{\frac{1-\varepsilon}{3}}|\Psi^{+}\rangle|10\rangle+\sqrt{\varepsilon}|\Psi^{-}\rangle|11\rangle$
(139)

Now for the rest of the circuit, a measurement-controlled system will be
applied based on the four orthogonal projective measurements. If we measured
$|00\rangle,|\Psi^{+}\rangle$ with probabilities $p_{1}$ and $p_{2}$, the
control system will do nothing and leave the principal system free to output
$\sigma_{1,2}$ guaranteed by Schmidt decomposition. Once the state
$|11\rangle$ is measured with probability $p_{4}$, a
$U_{\varepsilon}^{\dagger}$ will be applied to the main system in order to
erase the factor of $\varepsilon$ and produce $\sigma_{4}$. At last, if
$|\Psi^{-}\rangle$ is measured with probability $p_{3}$, a unitary
$W=(V_{\Psi^{-}}U^{\dagger}U_{\varepsilon}^{\dagger}\otimes I)CV^{\dagger}$
will be applied on both principal and the ancillary systems to reverse to the
initial state $|00\rangle|00\rangle$ and then generate a pure state
$\sigma_{3}=|\Psi^{-}\rangle\\!\langle\Psi^{-}|$.

##  Appendix H SDP for inverse operation method

Here is the SDP for finding the minimum overhead protocol of inverse operation
method.

$\displaystyle\min$ $\displaystyle\quad p$ (140) subject to
$\displaystyle\quad J_{\cal D}\geq 0$ (141)
$\displaystyle\quad\operatorname{Tr}_{B}[J_{{\cal D}_{AB}}]=pI_{A}$ (142)
$\displaystyle\quad J_{{\cal F}_{AC}}\equiv\operatorname{Tr}_{B}[(J_{{\cal
N}_{AB}}^{T_{B}}\otimes I_{C})(I_{A}\otimes J_{{\cal D}_{BC}})]$ (143)
$\displaystyle\quad{\cal F}_{AB}=J_{{\operatorname{id}}_{AB}}$ (144)

$J_{\cal D}$ and ${\cal N}$ are the Choi matrices for retriever ${\cal D}$ and
noisy ${\cal N}$ respectively. $J_{{\operatorname{id}}_{AB}}$ refers to the
Choi matrix for the identity channel. Eq. 141 and Eq. 142 corresponds to the
condition that the retriever ${\cal D}$ is CPTP. In Eq. 143, $J_{\cal F}$ is
the Choi matrix of the composed map ${\cal D}\circ{\cal N}$. Eq. 144
corresponds to the condition ${\cal D}={\cal N}^{-1}$.

##  Appendix I Proof for Lemma 6

$\displaystyle f_{\min}({\cal N},k)$ $\displaystyle=\min\Big{\\{}f\mid({\cal
N}^{\otimes k})^{\dagger}\circ\tilde{{\cal
C}}^{\dagger}(H)=\frac{1}{f}\left(H+tI\right),\tilde{{\cal
C}}\in\text{CPTP}\Big{\\}}$ (145) $\displaystyle g_{\min}({\cal N})$
$\displaystyle=\min\Big{\\{}\sum_{i}|c_{i}|\mid{\cal
N}^{-1}=\sum_{i}c_{i}{\cal C}_{i},{\cal C}_{i}\in\text{CPTP}\Big{\\}},$ (146)

Lemma 6 For arbitrary invertible quantum noisy channel ${\cal N}$, and moment
order $k$, we have

$f_{\min}({\cal N},k)\leq g_{\min}({\cal N},k).$ (147)

###### Proof.

In the task of extracting $k$-th moment from noisy states ${\cal N}(\rho)$, we
can apply the inverse operation with optimal decomposition on quantum systems
simultaneously to mitigate the error. The corresponding sampling overhead for
$k$-th moment will be

$g_{\min}({\cal N},k)=\min\Big{\\{}\sum_{i}|c_{i}|\mid({\cal N}^{-1})^{\otimes
k}=\sum_{i}c_{i}{\cal C}_{i},c_{i}\in\mathbb{R},{\cal
C}_{i}\in\text{CPTP}\Big{\\}}.$ (148)

Let’ assume the optimal decomposition with lowest sampling overhead is

$({\cal N}^{-1})^{\otimes k}=\sum_{i}c_{i}^{\prime}{\cal C}_{i}^{\prime}.$
(149)

Since $({\cal N}^{-1})^{\otimes k}\circ{\cal N}^{\otimes k}=I$, where $I$
refers to identity. we directly have ${\cal N}^{\otimes k\dagger}\circ({\cal
N}^{-1})^{\otimes k\dagger}(H)=H$, where $H$ is the moment observable. For
further step, we have ${\cal N}^{\otimes k\dagger}\circ({\cal
N}^{-1})^{\otimes k\dagger}(H+tI-tI)=H$, where $t$ is a real coefficient. Note
that ${\cal N}^{\dagger}$ is a CPUP map, and $({\cal N}^{-1})^{\dagger}$ is a
HPUP map, thus

${\cal N}^{\otimes k\dagger}\circ[({\cal N}^{-1})^{\otimes
k\dagger}(H)+tI]=H+tI.$ (150)

We can denote ${\cal C}^{\dagger}=[({\cal N}^{-1})^{\otimes k\dagger}(H)+tI]$
can be consider as a whole. If we take the coefficient $t$ greater than the
absolute value of the minimum eigenvalue of the hermitian matrix $({\cal
N}^{-1})^{\otimes k\dagger}(H)$, which is $t\geq|\min(\text{eig}({\cal
N}^{-1})^{\otimes k\dagger}(H))|$. Then, the map ${\cal C}^{\dagger}$ is a
CPUS map, and its corresponding adjoint map ${\cal C}$ is a CPTS map. We can
substitute the constructed map ${\cal C}$ back to Eq. (150), we have

$({\cal N}^{\otimes k})^{\dagger}\circ{\cal C}^{\dagger}(H)=H+tI$ (151)

Since any CPTS map ${\cal C}$ can be written in the form of ${\cal
C}=f\tilde{{\cal C}}$, where $f$ is the real coefficient, and $\tilde{{\cal
C}}$ is a CPTP map. We have

$({\cal N}^{\otimes k})^{\dagger}\circ f\tilde{{\cal C}}^{\dagger}(H)=H+tI,$
(152)

which has the same form as Eq. (3). Therefore, we have proved that the optimal
decomposition of Eq. (V) is one feasible solution of Eq. 3, and we get the
relation $f_{\min}({\cal N},k)\leq g_{\min}({\cal N},k)$ straightforwardly.
The proof is completed  $\sqcap$$\sqcup$

##  Appendix J DE retriever when $\rho$ is n-qubit state

###### Proposition 11

Given noised states ${\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes 2}$, and
error tolerance $\delta$ where $\rho$ is an $n$-qubit state, the second order
moment $\operatorname{Tr}[\rho^{2}]$ can be estimated by
$f\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes
2}]-t$ where $f=\frac{1}{(1-\epsilon)^{2}}$ and
$t=\frac{1-(1-\epsilon)^{2}}{2^{n}(1-\epsilon)^{2}}$. The term
$\operatorname{Tr}[H{\cal C}\circ{\cal N}^{\epsilon}_{\rm DE}(\rho)^{\otimes
2}]$ can be estimated by implementing quantum retriever ${\cal C}$ with
optimal sample compleixty ${\cal O}(1/\delta^{2}(1-\epsilon)^{4})$. The Choi
matrix of such a quantum retriever is,

$J_{\cal C}=\frac{1}{4^{n}}I^{\otimes n}I^{\otimes
n}+\frac{1}{4^{n}(4^{n}-1)}(\sum_{i\neq 0}P_{i}\otimes
P_{i})\otimes(\sum_{j\neq 0}P_{j}\otimes P_{j}).$

###### Proof.

For an arbitrary state $\rho\otimes\rho$, after applying the noise ${\cal N}$
onto it, we denote the state as $\rho^{\prime}={\cal N}(\rho\otimes\rho)$,
thus

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$ $\displaystyle={\cal
D}(\rho^{\prime})=\operatorname{Tr}_{A}[(\rho^{\prime T}\otimes I^{\otimes
n})J_{\cal C}]$ (153) $\displaystyle=\operatorname{Tr}_{A}\Big{[}(\rho^{\prime
T}\otimes I^{\otimes n})(\frac{1}{4^{n}}I^{\otimes n}I^{\otimes
n}+\frac{1}{4^{n}(4^{n}-1)}(\sum_{i\neq 0}P_{i}\otimes
P_{i})\otimes(\sum_{j\neq 0}P_{j}\otimes P_{j}))\Big{]}$ (154)
$\displaystyle=\frac{1}{4^{n}}\operatorname{Tr}[\rho^{\prime T}]I^{\otimes
n}+\frac{1}{4^{n}(4^{n}-1)}\operatorname{Tr}[\rho^{\prime T}(\sum_{i\neq
0}P_{i}\otimes P_{i})](\sum_{j\neq 0}P_{j}\otimes P_{j})$ (155)
$\displaystyle=\frac{1}{4^{n}}I^{\otimes
n}+\frac{1}{4^{n}(4^{n}-1)}\operatorname{Tr}[\rho^{\prime T}(\sum_{i\neq
0}P_{i}\otimes P_{i})](\sum_{j\neq 0}P_{j}\otimes P_{j})$ (156)

Note that the above equations utilized the face that transpose operation is
trace-preserving, i.e.,
$\operatorname{Tr}[\rho^{T}]=\operatorname{Tr}[\rho]=1$. Since the matrix
$\sum_{i\neq 0}P_{i}$ is symmetric, thus we have

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$
$\displaystyle=\frac{1}{4^{n}}I^{\otimes
n}+\frac{1}{4^{n}(4^{n}-1)}\operatorname{Tr}[\rho^{\prime T}(\sum_{i\neq
0}P_{i}\otimes P_{i})^{T}](\sum_{j\neq 0}P_{j}\otimes P_{j})$ (157)
$\displaystyle=\frac{1}{4^{n}}I^{\otimes
n}+\frac{1}{4^{n}(4^{n}-1)}\operatorname{Tr}[\rho^{\prime}(\sum_{i\neq
0}P_{i}\otimes P_{i})](\sum_{j\neq 0}P_{j}\otimes P_{j})$ (158)

The trace term can be calculated by substituting Eq. (80), we have

$\displaystyle\operatorname{Tr}[\rho^{\prime}(\sum_{i\neq 0}P_{i}\otimes
P_{i})]$
$\displaystyle=\operatorname{Tr}[((1-\epsilon)^{2}\rho\otimes\rho+\frac{\epsilon(1-\epsilon)}{2^{n}}I^{\otimes
n}\otimes\rho+\frac{(1-\epsilon)\epsilon}{2^{n}}\rho\otimes I^{\otimes
n}+\frac{\epsilon^{2}}{4^{n}}(I^{\otimes n}\otimes I^{\otimes n}))(\sum_{i\neq
0}P_{i}\otimes P_{i})]$ (159)
$\displaystyle=(1-\epsilon)^{2}\operatorname{Tr}[\rho\otimes\rho(\sum_{i\neq
0}P_{i}\otimes P_{i})\otimes(\sum_{j\neq 0}P_{j}\otimes P_{j})].$ (160)

In the second equation, all other terms are traceless, thus, only the first
term survives. Replace the equation back to Eq. (158), then

$\displaystyle{\cal C}\circ{\cal N}(\rho\otimes\rho)$
$\displaystyle=\frac{1}{4^{n}}I^{\otimes
n}+\frac{(1-\epsilon)^{2}}{4^{n}(4^{n}-1)}\operatorname{Tr}[\rho\otimes\rho(\sum_{i\neq
0}P_{i}\otimes P_{i})](\sum_{j\neq 0}P_{j}\otimes P_{j}).$ (161)

The information $\operatorname{Tr}[\rho^{2}]$ is estimated from
$\operatorname{Tr}[H(\rho\otimes\rho)]$, where $H=\frac{1}{2^{n}}\sum_{i\neq
0}P_{i}\otimes P_{i}$ is cyclic permutation operator. It is easy to check that

$\displaystyle\operatorname{Tr}[H{\cal C}\circ{\cal N}(\rho\otimes\rho)]$
$\displaystyle=\frac{1}{4^{n}}\operatorname{Tr}[H*I^{\otimes
n}]+\frac{(1-\epsilon)^{2}}{4^{n}(4^{n}-1)}\operatorname{Tr}[\rho\otimes\rho(\sum_{i\neq
0}P_{i}\otimes P_{i})]\operatorname{Tr}[H*(\sum_{j\neq 0}P_{j}\otimes
P_{j})].$ (162)

We can quickly get $\operatorname{Tr}[H]=2^{n}$ and
$\operatorname{Tr}[H*(\sum_{j\neq 0}P_{j}\otimes P_{j})]=2^{n}(4^{n}-1)$. Then

$\displaystyle f\operatorname{Tr}[H{\cal C}\circ{\cal N}(\rho\otimes\rho)]$
$\displaystyle=\frac{1}{(1-\epsilon)^{2}}\Big{[}\frac{1}{2^{n}}+\frac{(1-\epsilon)^{2}}{2^{n}}\operatorname{Tr}[\rho\otimes\rho(\sum_{i\neq
0}P_{i}\otimes P_{i})]\Big{]}$ (163)
$\displaystyle=\frac{1}{2^{n}(1-\epsilon)^{2}}+\frac{1}{2^{n}}\operatorname{Tr}[\rho\otimes\rho(\sum_{i\neq
0}P_{i}\otimes P_{i})]$ (164)
$\displaystyle=\frac{1}{2^{n}(1-\epsilon)^{2}}-\frac{1}{2^{n}}+\frac{1}{2^{n}}+\frac{1}{2^{n}}\operatorname{Tr}[\rho\otimes\rho(\sum_{i\neq
0}P_{i}\otimes P_{i})]$ (165)
$\displaystyle=\frac{1}{2^{n}(1-\epsilon)^{2}}-\frac{(1-\epsilon)^{2}}{2^{n}(1-\epsilon)^{2}}+\frac{1}{2^{n}}\operatorname{Tr}[\rho\otimes\rho(\sum_{i}P_{i}\otimes
P_{i})]$ (166)
$\displaystyle=\frac{1-(1-\epsilon)^{2}}{2^{n}(1-\epsilon)^{2}}+\operatorname{Tr}[\rho\otimes\rho
H]$ (167) $\displaystyle=t+\operatorname{Tr}[\rho^{2}].$ (168)

where $t=\frac{1-(1-\epsilon)^{2}}{2^{n}(1-\epsilon)^{2}}$ is the shift
coefficient. The desired high-order moment $\operatorname{Tr}[\rho^{2}]$ is
directly obtained by estimated value $\operatorname{Tr}[H{\cal C}\circ{\cal
N}(\rho\otimes\rho)]$ substrate a constant $t$.  $\sqcap$$\sqcup$