3}(|b\rangle\langle b|^{\otimes
3})\left(\mathcal{M}^{-1}(O_{i})\otimes\mathcal{M}^{-1}(O_{i}^{\dagger})\otimes\rho\right)\right]-|\mathrm{tr}(O_{i}\rho)|^{2},$
(231)

using the fact that $\mathcal{M}$ is self-adjoint (with respect to the
Hilbert-Schmidt inner product) in the last line, so
$\mathrm{tr}(\mathcal{M}^{-1}(A)B)=\mathrm{tr}(A\mathcal{M}^{-1}(B))$ for any
operators $A$ and $B$. This variance can be upper-bounded by the first term,
maximised over all states $\rho$; this gives the square of what is often
referred to as the “shadow norm” in the literature. A useful observation is
that since $\mathrm{tr}(\hat{\rho})=1$, the variance for any observable
$\hat{O}_{i}$ is the same as that for $\hat{O}_{i}+cI$ for any
$c\in\mathbb{C}$, so we can also also write

$\displaystyle\mathrm{Var}[\hat{o}_{i}]=\mathrm{tr}\left[\sum_{b\in\mathcal{B}}\mathop{{}\mathbb{E}}_{U\in\mathcal{D}}\mathcal{U}^{\otimes
3}(|b\rangle\langle b|^{\otimes
3})\left(\mathcal{M}^{-1}(O_{i}+cI)\otimes\mathcal{M}^{-1}(O_{i}^{\dagger}+cI)\otimes\rho\right)\right]-|\mathrm{tr}((O_{i}+cI)\rho)|^{2}.$
(232)

If median-of-means estimators are used, it follows straightforwardly from
Chebyshev’s and Hoeffding’s inequalities that

$N_{\mathrm{sample}}=\mathcal{O}\left(\frac{\log(M/\delta)}{\varepsilon^{2}}\max\limits_{i\in[M]}\mathrm{Var}[\hat{o}_{i}]\right)$
(233)

classical shadows samples (and hence single-copy measurements of $\rho$) are
sufficient to ensure that with probability at least $1-\delta$, every
$\mathrm{tr}(O_{i}\rho)$ is estimated to within additive error $\varepsilon$.

###  H.2 Preliminaries and notation

####  H.2.1 Generalized Clifford group

The generalized (single-qudit) Clifford group in dimension $d$, which we will
denote by $\mathrm{Cl}_{d}$, is the normaliser of the Pauli group in the
$d$-dimensional unitary group.222Strictly speaking, we should quotient out
phases so that it makes sense to write the twirl channel in Eq. (238) as an
average over a finite group, but the results in this appendix hold even if
phases are included. Hence, a Clifford $U\in\mathrm{Cl}_{d}$ satisfies the
property that for any $(q,p)\in\mathbb{Z}_{d}^{2}$,
$U^{\dagger}D_{q,p}U=\alpha D_{q^{\prime},p^{\prime}}$ for some
$(q^{\prime},p^{\prime})\in\mathbb{Z}_{d}^{2}$ and some phase $\alpha$. Thus,
each Clifford can be associated with a map
$\mathbb{Z}_{d}^{2}\to\mathbb{Z}_{d}^{2}$, and since its action under
conjugation is completely specified by how $X=D_{1,0}$ and $Z=D_{0,1}$ are
transformed, this map is linear, and can be written as a $2\times 2$ matrix
$C\in\mathbb{Z}_{d}^{2\times 2}$ so that

$U^{\dagger}D_{q,p}U=\alpha D_{q^{\prime},p^{\prime}}\text{ for some phase
{\hbox{\alpha}} }\quad\Leftrightarrow\quad C\begin{pmatrix}q\\\
p\end{pmatrix}=\begin{pmatrix}q^{\prime}\\\ p^{\prime}\end{pmatrix}$ (234)

for all $q,p\in\mathbb{Z}_{d}$. It follows from the fact that conjugation
preserves commutation relations that this matrix also has to be symplectic. It
turns out that conversely, for any symplectic $2\times 2$ matrix with entries
in $\mathbb{Z}_{d}$, there exists a corresponding Clifford in the sense of Eq.
(234) [50].

######  Fact 1 ([50]).

For any $C\in\mathbb{Z}_{d}^{2\times 2}$, there exists $U\in\mathrm{Cl}_{d}$
satisfying Eq. (234) if and only if $C$ is symplectic.

The following basic fact about symplectic $2\times 2$ matrices will also come
in useful.

######  Fact 2.

A $2\times 2$ matrix $C$ is symplectic if and only if $\det(C)=1$.

####  H.2.2 Liouville representation

In some parts of this appendix (especially when we analyse the twirl channels
of the generalized Clifford group in the following subsection), we will use
Liouville representation, which is often useful for specifying quantum
channels. In this representation, operators are notated using “double” kets,
and a “double” bracket is used to represent the Hilbert-Schmidt inner product.
By convention, all double kets are normalised with respect to the Hilbert-
Schmidt norm, unless otherwise noted. Thus, for any nonzero operators $A,B$,

$|A\rrangle\equiv\frac{1}{\sqrt{\mathrm{tr}(A^{\dagger}A)}}A,\qquad\llangle
A|B\rrangle\equiv\frac{\mathrm{tr}(A^{\dagger}B)}{\sqrt{\mathrm{tr}(A^{\dagger}A)\mathrm{tr}(B^{\dagger}B)}},$
(235)

and we set $|0\rangle\equiv 0$ for the zero operator. In particular, since the
generalized Pauli operators are unitary (Eq. (13)) and Hilbert-Schmidt
orthogonal (Proposition 8), we have

$|D_{q,p}\rrangle\equiv\frac{1}{\sqrt{d}}D_{q,p},\qquad\llangle
D_{q,p}|D_{q^{\prime},p^{\prime}}\rrangle=\delta_{q,q^{\prime}}\delta_{p,p^{\prime}}.$
(236)

As with usual (state) kets, we will freely write e.g., $|A\rrangle|B\rrangle$
in place of $|A\rrangle\otimes|B\rrangle$.

A superoperator $\mathcal{E}$ acting on an operator $A$ is represented by
placing the superoperator to the left of the operator’s double ket:

$\mathcal{E}|A\rrangle\equiv\frac{1}{\sqrt{\mathrm{tr}(A^{\dagger}A)}}\mathcal{E}(A).$

We will also write $\mathcal{E}\mathcal{E}^{\prime}$ in place of
$\mathcal{E}\circ\mathcal{E}^{\prime}$.

Since the generalized Pauli operators form a Hilbert-Schmidt orthogonal basis
for the space of operators acting on $d$-dimensional Hilbert space
$\mathcal{H}_{d}$, we have the following resolution of the identity
superoperator $\mathcal{I}$:

$\mathcal{I}=\sum_{p,q\in\mathbb{Z}_{d}}|D_{p,q}\rrangle\llangle D_{p,q}|.$
(237)

###  H.3 Uniform distribution over generalized Cliffords

In the rest of this appendix, we consider the classical shadows associated
with the uniform distribution over $\mathrm{Cl}_{d}$, the generalized (single-
qudit) Cliffords in prime dimension $d$.

As can be seen from Eqs. (225) and (231), the measurement channel
$\mathcal{M}$ in the classical shadows protocol depends on the chosen unitary
distribution $\mathcal{D}$ only through its $2$-fold twirl
$\mathop{{}\mathbb{E}}_{U\sim\mathcal{D}}\mathcal{U}^{\otimes 2}$, while the
variances of the resulting estimates depend on $\mathcal{D}$ through its
$3$-fold twirl $\mathop{{}\mathbb{E}}_{U\sim\mathcal{D}}\mathcal{U}^{\otimes
3}$. Hence, in this subsection, we evaluate the $2$\- and $3$-fold twirl
channels for the uniform distribution over generalized Cliffords, arriving at
the theorem below. This then allow us to obtain explicit expressions for the
measurement channel $\mathcal{M}$ and the variance for any observable.

For $k\in\mathbb{Z}_{>0}$, we use $\mathcal{E}_{d}^{(k)}$ to denote the
$k$-fold twirl channel corresponding to the uniform distribution over
$\mathrm{Cl}_{d}$:

$\mathcal{E}_{d}^{(k)}\coloneqq\frac{1}{|\mathrm{Cl}_{d}|}\sum_{U\in\mathrm{Cl}_{d}}\mathcal{U}^{\otimes
k},$ (238)

where, as in subsection H.1, $\mathcal{U}(\,\cdot\,)=U^{\dagger}(\,\cdot\,)U$.

######  Theorem 25 (First three moments of uniform distribution over
generalized Cliffords).

Let $\mathcal{E}_{d}^{(k)}$ be defined as in Eq. (238). Then, for any prime
$d$, we have

  1. (i)

$\mathcal{E}_{d}^{(1)}=|I\rrangle\llangle I|,$

  2. (ii)

$\mathcal{E}_{d}^{(2)}=|I\rrangle\llangle I|^{\otimes
2}+|\Phi_{0}\rrangle\llangle\Phi_{0}|$,

  3. (iii)

$\mathcal{E}_{d}^{(3)}=|I\rrangle\llangle I|^{\otimes
3}+|\Phi_{1}\rrangle\llangle\Phi_{1}|+|\Phi_{2}\rrangle\llangle\Phi_{2}|+|\Phi_{3}\rrangle\llangle\Phi_{3}|+\sum\limits_{k=1}^{d-2}|\Psi_{k}\rrangle\llangle\Psi_{k}|+\sum\limits_{l=1}^{d-1}|\Upsilon_{l}\rrangle\llangle\Upsilon_{l}|$,

where

$|\Phi_{0}\rrangle\coloneqq\frac{1}{\sqrt{d^{2}-1}}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle$ (239)

and

$\displaystyle|\Phi_{1}\rrangle\coloneqq\frac{1}{\sqrt{d^{2}-1}}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|I\rrangle|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle,$
(240)
$\displaystyle|\Phi_{2}\rrangle\coloneqq\frac{1}{\sqrt{d^{2}-1}}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|I\rrangle|D_{q,p}^{\dagger}\rrangle,$
(241)
$\displaystyle|\Phi_{3}\rrangle\coloneqq\frac{1}{\sqrt{d^{2}-1}}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle|I\rrangle,$
(242)
$\displaystyle|\Psi_{k}\rangle\coloneqq\frac{1}{\sqrt{d^{2}-1}}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|D_{kq,kp}\rrangle|D_{q,p}^{\dagger}D_{kq,kp}^{\dagger}\rrangle$
(243)
$\displaystyle|\Upsilon_{l}\rrangle\coloneqq\frac{1}{\sqrt{d(d^{2}-1)}}\sum_{\begin{subarray}{c}q_{1},p_{1},q_{2},p_{2}\in\mathbb{Z}_{d}\\\
p_{1}q_{2}-q_{1}p_{2}=l\end{subarray}}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}\rrangle.$
(244)

The strategy we use to prove Theorem 25 is very similar in spirit to that used
to determine the $2$\- and $3$-fold twirl channels of the matchgate group in
Ref. [31]. In broad strokes, we start from the simple observation that each
twirl channel is an orthogonal projector, and use the invariance of the twirl
channel under composition with suitably chosen group elements to determine its
image. The difference in the details of the proof of Theorem 25 differ from
those in Ref. [31] stem from the difference in the structure of the adjoint
representation of the Clifford group and that of the matchgate group, with the
former being somewhat easier to handle since it is a discrete group.

The following facts follow straightforwardly from the fact that (for any
$k\in\mathbb{Z}_{>0}$) the map $U\mapsto\mathcal{U}^{\otimes k}$ is a
representation of the generalized Clifford group; see e.g., [51] Proposition
2.8.

######  Fact 3.

For any $k\in\mathbb{Z}_{>0}$, $\mathcal{E}_{d}^{(k)}$ is the orthogonal
projector onto the subspace $\\{A\in\mathcal{L}(\mathcal{H}_{d})^{\otimes
k}:\mathcal{U}^{\otimes k}(A)=A\enspace\forall\,U\in\mathrm{Cl}_{d}\\}$.

######  Fact 4.

For any $k\in\mathbb{Z}_{>0}$ and $U\in\mathrm{Cl}_{d}$,

$\mathcal{U}^{\otimes
k}\circ\mathcal{E}_{d}^{(k)}=\mathcal{E}_{d}^{(k)}=\mathcal{E}_{d}^{(k)}\circ\mathcal{U}^{\otimes
k}.$

We will also make use of the following simple fact.

######  Fact 5.

Let $d$ be a prime number. For any $q,p\in\mathbb{Z}_{d}$ with
$(q,p)\neq(0,0)$, there exist $q^{\prime},p^{\prime}\in\mathbb{Z}_{d}$ such
that

$p^{\prime}q-q^{\prime}p=1.$

###### Proof.

If $q\neq 0$, choose any $q^{\prime}$ and set
$p^{\prime}=(q^{\prime}p+1)q^{-1}$. Otherwise, we must have $p\neq 0$, in
which case we choose any $p^{\prime}$ and set
$q^{\prime}=(p^{\prime}q-1)p^{-1}$. ∎

####  H.3.1 1-fold twirl

We start by evaluating the $1$-fold twirl $\mathcal{E}_{d}^{(1)}$, since a
similar idea will be used in obtaining the $2$\- and $3$-fold twirls.

###### Proof of Theorem 25(i).

We can rewrite Eq. (16) as

$D_{q^{\prime},p^{\prime}}D_{q,p}D_{q^{\prime},p^{\prime}}^{\dagger}=\omega^{p^{\prime}q-q^{\prime}p}D_{q,p},$
(245)

with $\omega\coloneqq e^{2\pi i/d}$. By Fact 5, there exist
$q^{\prime},p^{\prime}$ such that $p^{\prime}q-q^{\prime}p=1$ provided that
$(q,p)\neq(0,0)$. Thus, for any $(q,p)\neq(0,0)$, there exists a Pauli
$U\in\mathrm{Cl}_{d}$ such that

$\mathcal{U}|D_{q,p}\rrangle=\omega|D_{q,p}\rrangle.$

Hence,

$\displaystyle\mathcal{E}_{d}^{(1)}|D_{q,p}\rrangle$
$\displaystyle=\frac{1}{d}\sum_{j=0}^{d-1}\mathcal{E}_{d}^{(1)}\mathcal{U}^{j}|D_{q,p}\rrangle$
$\displaystyle=\frac{1}{d}\mathcal{E}_{d}^{(1)}\left(\sum_{j=0}^{d-1}\omega^{j}\right)|D_{q,p}\rrangle$
$\displaystyle=0$

for any $(q,p)\neq(0,0)$, where we use Fact 4 in the first line, while for
$(q,p)=(0,0)$, we have
$\mathcal{E}_{d}^{(1)}|D_{0,0}\rrangle=\mathcal{E}_{d}^{(1)}|I\rrangle=|I\rrangle$.
Since $\\{D_{q,p}\\}_{q,p\in\mathbb{Z}_{d}}$ forms a basis for the space of
operators, it follows that

$\mathcal{E}_{1}=|I\rrangle\llangle I|.$

∎

####  H.3.2 2-fold twirl

Now, we evaluate the $2$-fold twirl $\mathcal{E}_{d}^{(2)}$, which will allow
us to determine the measurement channel (Eq. (225)) for the generalized
Clifford shadows. We start with a lemma that precludes certain basis states
(in the Pauli basis) from being in the image of $\mathcal{E}_{d}^{(2)}$.

######  Lemma 26.

Let $d$ be a prime number. For $q_{1},p_{1},q_{2},p_{2}\in\mathbb{Z}_{d}$,
$\mathcal{E}_{d}^{(2)}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle\neq 0$
only if $(q_{1}+q_{2},p_{1}+p_{2})=(0,0)$.

###### Proof.

The proof is very similar to that of Theorem 25(i) above. From Eq. (245),

$D_{q^{\prime},p^{\prime}}^{\otimes 2}(D_{q_{1},p_{1}}\otimes
D_{q_{2},p_{2}})(D_{q^{\prime},p^{\prime}}^{\otimes
2})^{\dagger}=\omega^{p^{\prime}(q_{1}+q_{2})-q^{\prime}(p_{1}+p_{2})}D_{q_{1},p_{1}}\otimes
D_{q_{2},p_{2}},$

so it follows from Fact 5 that for any
$q_{1},p_{1},q_{2},p_{2}\in\mathbb{Z}_{d}$ with
$(q_{1}+q_{2},p_{1}+p_{2})\neq(0,0)$, there exists a Pauli
$U\in\mathrm{Cl}_{d}$ such that

$\mathcal{U}^{\otimes
2}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle=\omega|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle.$

Hence,

$\displaystyle\mathcal{E}_{d}^{(2)}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle$
$\displaystyle=\frac{1}{d}\sum_{j=0}^{d-1}\mathcal{E}_{d}^{(2)}(\mathcal{U}^{j})^{\otimes
2}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle$
$\displaystyle=\frac{1}{d}\mathcal{E}_{d}^{(2)}\left(\sum_{j=0}^{d-1}\omega^{j}\right)|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle$
$\displaystyle=0.$

∎

Next, we show that all $D_{q,p}$ with $(q,p)\neq(0,0)$ are in the same orbit
under the action of $\mathrm{Cl}_{d}$ by conjugation (if we ignore global
phases).

######  Lemma 27.

Let $d$ be a prime number. For any
$q,p,q^{\prime},p^{\prime}\in\mathbb{Z}_{d}$ with
$(q,p),(q^{\prime},p^{\prime})\neq(0,0)$, there exists $U\in\mathrm{Cl}_{d}$
such that

$\mathcal{U}|D_{q,p}\rrangle=\alpha|D_{q^{\prime},p^{\prime}}\rrangle$

for some phase $\alpha$.

###### Proof.

By Fact 1 and Eq. (234), such a Clifford exists provided that there exists a
symplectic matrix $\begin{pmatrix}a&b\\\ c&d\end{pmatrix}$ such that

$\begin{pmatrix}w&x\\\ y&z\end{pmatrix}\begin{pmatrix}q\\\
p\end{pmatrix}=\begin{pmatrix}q^{\prime}\\\ p^{\prime}\end{pmatrix}.$

If $q\neq 0$, we set $w=(q^{\prime}-xp)q^{-1}$ and $y=(p^{\prime}-zp)q^{-1}$.
Then, $\det\begin{pmatrix}w&x\\\
y&z\end{pmatrix}=(q^{\prime}z-xp^{\prime})q^{-1}$. By Fact 5, since
$(q^{\prime},p^{\prime})\neq(0,0)$, we can choose $z,x$ such that
$(q^{\prime}z-xp^{\prime})=q$, so that the determinant is $1$—i.e., the matrix
is symplectic by Fact 2. Otherwise, we must have $p\neq 0$, in which case we
set $x=(q^{\prime}-wq)p^{-1}$ and $z=(p^{\prime}-yq)p^{-1}$, and choose $w,y$
such that $(wp^{\prime}-q^{\prime}y)=p$ via Fact 5. ∎

###### Proof of Theorem 25(ii).

Inserting resolutions of the identity (Eq. (237)) and using Lemma 26, we have

$\displaystyle\mathcal{E}_{d}^{(2)}$
$\displaystyle=\sum_{q,p,q^{\prime},p^{\prime}\in\mathbb{Z}_{d}}|D_{q,p}\rrangle|D_{-q,-p}\rrangle\llangle
D_{q,p}|\llangle
D_{-q,-p}|\mathcal{E}_{d}^{(2)}|D_{q^{\prime},p^{\prime}}\rrangle|D_{-q^{\prime},-p^{\prime}}\rrangle\llangle
D_{q^{\prime},p^{\prime}}|\llangle D_{-q^{\prime},-p^{\prime}}|$ (246)
$\displaystyle=\sum_{q,p,q^{\prime},p^{\prime}\in\mathbb{Z}_{d}}|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle\llangle
D_{q,p}|\llangle
D_{q,p}^{\dagger}|\mathcal{E}_{d}^{(2)}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle\llangle
D_{q^{\prime},p^{\prime}}|\llangle D_{q^{\prime},p^{\prime}}^{\dagger}|$ (247)
$\displaystyle=|I\rrangle|I\rrangle\llangle I|\llangle
I|+\sum_{\begin{subarray}{c}q,p,q^{\prime},p^{\prime}\in\mathbb{Z}_{d}\\\
(q,p),(q^{\prime},p^{\prime})\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle\llangle
D_{q,p}|\llangle
D_{q,p}^{\dagger}|\mathcal{E}_{d}^{(2)}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle\llangle
D_{q^{\prime},p^{\prime}}|\llangle D_{q^{\prime},p^{\prime}}^{\dagger}|$ (248)

(note that since $\mathcal{E}_{d}^{(2)\dagger}=\mathcal{E}_{d}^{(2)}$ by Fact
3, Lemma 26 also implies that $\llangle D_{q_{1},p_{1}}|\llangle
D_{q_{2},p_{2}}|\mathcal{E}_{d}^{(2)}=0$ if
$(q_{2},p_{2})\neq(-q_{1},-p_{1})$.) Here, the second line uses Eq. (13) and
the third line follows from the obvious fact that
$\mathcal{E}_{d}^{(2)}|D_{0,0}\rrangle=|D_{0,0}\rrangle=|I\rrangle$.

Now, we use Lemma 27 to show that $\llangle D_{q,p}|\llangle
D_{q,p}^{\dagger}|\mathcal{E}_{d}^{(2)}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle$
is the same for all $(q,p),(q^{\prime},p^{\prime})\neq(0,0)$, i.e.,

$\llangle D_{q,p}|\llangle
D_{q,p}^{\dagger}|\mathcal{E}_{d}^{(2)}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle=c$
(249)

for some constant $c$. To see this, note from Lemma 27 that for any
$(q^{\prime\prime},p^{\prime\prime})\neq(0,0)$, there exists
$U\in\mathrm{Cl}_{d}$ such that
$\mathcal{U}|D_{q^{\prime},p^{\prime}}\rrangle=\alpha|D_{q^{\prime\prime},p^{\prime\prime}}\rrangle$
for some phase $\alpha$, so

$\mathcal{U}^{\otimes
2}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle=\alpha|D_{q^{\prime\prime},p^{\prime\prime}}\rrangle\alpha^{*}|D_{q^{\prime\prime},p^{\prime\prime}}^{\dagger}\rrangle=|D_{q^{\prime\prime},p^{\prime\prime}}\rrangle|D_{q^{\prime\prime},p^{\prime\prime}}^{\dagger}\rrangle$
(250)

Then, use Fact 4 to write

$\displaystyle\mathcal{E}_{d}^{(2)}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle$
$\displaystyle=\mathcal{E}_{d}^{(2)}\mathcal{U}^{\otimes
2}|D_{q^{\prime},p^{\prime}}\rrangle|D_{q^{\prime},p^{\prime}}^{\dagger}\rrangle=\mathcal{E}_{d}^{(2)}|D_{q^{\prime\prime},p^{\prime\prime}}\rrangle|D_{q^{\prime\prime},p^{\prime\prime}}^{\dagger}\rrangle.$

Inserting Eq. (249) into Eq. (248), we arrive at

$\displaystyle\mathcal{E}_{d}^{(2)}$ $\displaystyle=|I\rrangle\llangle
I|^{\otimes
2}+c\sum_{\begin{subarray}{c}q,p,q^{\prime},p^{\prime}\in\mathbb{Z}_{d}\\\
(q,p),(q^{\prime},p^{\prime})\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle\llangle
D_{q^{\prime},p^{\prime}}|\llangle D_{q^{\prime},p^{\prime}}^{\dagger}|$
$\displaystyle=|I\rrangle\llangle I|^{\otimes
2}+c|\Phi_{0}\rrangle\llangle\Phi_{0}|,$

where $|\Phi_{0}\rrangle$ is defined as in Eq. (239). Since
$\mathcal{E}_{d}^{(2)}$ is a projector (Fact 3), $c$ must be equal to $0$ or
$1$. We see from Eq. (250) that for any $U\in\mathrm{Cl}_{d}$,
$\mathcal{U}^{\otimes 2}$ permutes the basis states of the form
$|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle$ for $(q,p)\neq(0,0)$, without
incurring any phases. Since $|\Phi_{0}\rrangle$ is the equal superposition
over these basis states, it follows that $\mathcal{U}^{\otimes
2}|\Phi_{0}\rrangle=|\Phi_{0}\rrangle$, so
$\mathcal{E}_{d}^{(2)}|\Phi_{0}\rrangle=|\Phi_{0}\rrangle$ and $c$ must be
$1$. ∎

####  H.3.3 3-fold twirl

Finally, we evaluate the $3$-fold twirl $\mathcal{E}_{d}^{(3)}$, which will
allow us to analyse the variance (Eq. (231)) of estimates obtained using
generalized Clifford shadows. The high-level strategy for evaluating
$\mathcal{E}_{d}^{(3)}$ is is essentially the same as that for evaluating the
$2$-fold twirl channel $\mathcal{E}_{d}^{(2)}$. We start with the following
lemma, which precludes certain basis states from being in the image of
$\mathcal{E}_{d}^{(3)}$. This lemma is the straightforward extension of Lemma
26, which was used in evaluating $\mathcal{E}_{d}^{(2)}$. (In fact, there is a
natural generalization of these lemmas to $\mathcal{E}_{d}^{(k)}$ for any
$k\in\mathbb{Z}_{>0}$.)

######  Lemma 28.

Let $d$ be a prime number. For
$q_{1},p_{1},q_{2},p_{2},q_{3},p_{3}\in\mathbb{Z}_{d}$,
$\mathcal{E}_{d}^{(3)}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{3},p_{3}}\rrangle\neq
0$ only if $(q_{1}+q_{2}+q_{3},p_{1}+p_{2}+p_{3})=(0,0)$.

###### Proof.

The proof is essentially the same as that of Lemma 26, except we consider
three tensor components instead of two. ∎

Lemma 28 shows that we need only consider the action of
$\mathcal{E}_{d}^{(3)}$ on basis states of the form
$|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{-q_{1}-q_{2},-p_{1}-p_{2}}\rrangle$.
For these states, the third tensor component is fully determined by the first
two. For the purpose of the next lemma, we will say that two pairs of Paulis
$(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$ and
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ are in
the “orbit”333Strictly speaking, these are not orbits in the group theoretic
sense, as we allow phases in Eq. (251). However, this definition will suffice
for our purposes, because as can be seen from Eq. (253), we will eventually
choose the particular form of the basis states for
$\mathcal{L}(\mathcal{H}_{d})^{\otimes 3}$ in such a way that any phases from
Eq. (251) are cancelled out. if there exists $U\in\mathrm{Cl}_{d}$ such that
$\mathcal{U}^{\otimes
2}|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle=\alpha|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle$
for some phase $\alpha$, or equivalently,

$\mathcal{U}|D_{q_{1},p_{1}}\rrangle=\alpha_{1}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle\quad\text{
and
}\quad\mathcal{U}|D_{q_{2},p_{2}}\rrangle=\alpha_{2}|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle$
(251)

for some phases $\alpha_{1}$ and $\alpha_{2}$. The next lemma is essentially
the analogue of Lemma 27, but for pairs of Paulis.

######  Lemma 29.

The set of pairs of Paulis
$\\{(D_{q_{1},p_{1}},D_{q_{2},p_{2}}):q_{1},p_{1},q_{2},p_{2}\in\mathbb{Z}_{d}\\}$
is partitioned into the following orbits (in the sense of Eq. (251)):

  1. a)

$\\{(D_{0,0},D_{0,0})\\}$ ($1$ element)

  2. b)

$\\{(D_{0,0},D_{q,p}):q,p\in\mathbb{Z}_{d},(q,p)\neq(0,0)\\}$ ($d^{2}-1$
elements)

  3. c)

$\\{(D_{q,p},D_{0,0}):q,p\in\mathbb{Z}_{d},(q,p)\neq(0,0)\\}$ ($d^{2}-1$
elements),

  4. d)

$\\{(D_{q,p},D_{kq,kp}):q,p\in\mathbb{Z}_{d},(q,p)\neq(0,0)\\}$ for
$k\in\\{1,\dots,d-1\\}$ ($d^{2}-1$ elements for each $k$)

  5. e)

$\\{(D_{q_{1},p_{1}},D_{q_{2},p_{2}}):q_{1},p_{1},q_{2},p_{2}\in\mathbb{Z}_{d},p_{1}q_{2}-q_{1}p_{2}=l\\}$
for $l\in\\{1,\dots,d-1\\}$ ($d(d^{2}-1)$ elements for each $l$).

###### Proof.

First, we show that these subsets indeed partition the set of all pairs of
Paulis. The subsets in a), b) and c) are clearly mutually disjoint. The
remaining pairs are of the form $(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$ with
$(q_{1},p_{1}),(q_{2},p_{2})\neq 0$. The pairs with $p_{1}q_{2}-q_{1}p_{2}\neq
0$ can be partitioned into the subsets in e). This leaves pairs
$(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$ with $(q_{1},p_{1}),(q_{2},p_{2})\neq 0$
and $p_{1}q_{2}-q_{1}p_{2}=0$. Since $(q_{1},p_{1}),(q_{2},p_{2})\neq 0$, we
must have that $q_{1},q_{2}\neq 0$ and/or $p_{1},p_{2}\neq 0$. If
$q_{1},q_{2}\neq 0$, then we can write $q_{2}=(q_{1}^{-1}q_{2})q_{1}$, and
from $p_{1}q_{2}-q_{1}p_{2}=0$, we have $p_{2}=(q_{1}^{-1}q_{2})p_{1}$. If
$p_{1},p_{2}\neq 0$, then $q_{2}=(p_{1}^{-1}p_{2})q_{1}$ and
$p_{2}=(p_{1}^{-1}p_{2})p_{1}$. In both cases, $q_{2}=kq_{1}$ and
$p_{2}=kp_{1}$ for some $k\in\\{1,\dots,d-1\\}$, which coincides with the
subsets in d). It remains to check that these subsets are indeed orbits.

  1. a)

Since $\mathcal{U}|D_{0,0}\rrangle=|D_{0,0}\rrangle$ (recall $D_{0,0}=I$) for
any $U\in\mathrm{Cl}_{d}$, it is clear that $\\{(D_{0,0},D_{0,0})\\}$ is an
orbit.

  2. b)

Since $\mathcal{U}|D_{0,0}\rrangle=|D_{0,0}\rrangle$ for any
$U\in\mathrm{Cl}_{d}$,
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ is in
the same orbit as $(D_{0,0},D_{q,p})$ if and only if
$(q_{1}^{\prime},p_{1}^{\prime})=(0,0)$ and there exists $U\in\mathrm{Cl}_{d}$
such that
$\mathcal{U}|D_{q,p}\rrangle\propto|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle$.
By Lemma 27, for $(q,p)\neq(0,0)$, the latter is true whenever
$(q_{2}^{\prime},p_{2}^{\prime})\neq 0$. Hence,
$\\{(D_{0,0},D_{q,p}):q,p\in\mathbb{Z}_{d},(q,p)\neq(0,0)\\}$ is an orbit.

  3. c)

Same proof as for b).

  4. d)

For $(q,p)\neq(0,0)$ and $k\neq 0$,
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ is in
the same orbit as $(D_{q,p},D_{kq,kp})$ if and only if there exists
$U\in\mathrm{Cl}_{d}$ such that
$\mathcal{U}|D_{q,p}\rrangle=|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle$ and
$\mathcal{U}|D_{kq,kp}\rrangle=|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle$. By
Lemma 27, there exists $U\in\mathrm{Cl}_{d}$ such that
$\mathcal{U}|D_{q,p}\rrangle=|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle$
whenever $(q_{1}^{\prime},p_{1}^{\prime})\neq(0,0)$. Then, since
$D_{kq,kp}=D_{q,p}^{k}$, we have $\mathcal{U}(D_{kq,kp})\propto
U^{\dagger}(D_{q,p}^{k})U=D_{q_{1}^{\prime},p_{1}^{\prime}}^{k}=D_{kq_{1}^{\prime},kp_{1}^{\prime}}$.
Thus, all elements in the same orbit as $(D_{q,p},D_{kq,kp})$ are of the form
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{kq_{1}^{\prime},kp_{2}^{\prime}})$.

  5. e)

Finally, we show that for $p_{1}q_{2}-q_{1}p_{2}\neq 0$,
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ is in
the same orbit as $(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$—i.e., there exists
$U\in\mathrm{Cl}_{d}$ satisfying Eq. (251)—if and only if
$p_{1}^{\prime}q_{2}^{\prime}-q_{1}^{\prime}p_{2}^{\prime}=p_{1}q_{2}-q_{1}p_{2}$.
It follows from Facts 1 and 2 that there exists such a Clifford $U$ if and
only if there exists a $2\times 2$ matrix $C$ such that

$C\begin{pmatrix}q_{1}&q_{2}\\\
p_{1}&p_{2}\end{pmatrix}=\begin{pmatrix}q_{1}^{\prime}&q_{2}^{\prime}\\\
p_{1}^{\prime}&p_{2}^{\prime}\end{pmatrix}$ (252)

and $\det(C)=1$. Since $\det\begin{pmatrix}q_{1}&q_{2}\\\
p_{1}&p_{2}\end{pmatrix}=q_{1}p_{2}-p_{1}q_{2}\neq 0$ by assumption,
$\begin{pmatrix}q_{1}&q_{2}\\\ p_{1}&p_{2}\end{pmatrix}$ is invertible and Eq.
(252) is equivalent to $C=\begin{pmatrix}q_{1}^{\prime}&q_{2}^{\prime}\\\
p_{1}^{\prime}&p_{2}^{\prime}\end{pmatrix}\begin{pmatrix}q_{1}&q_{2}\\\
p_{1}&p_{2}\end{pmatrix}^{-1}$. Then,

$\det(C)=\det\begin{pmatrix}q_{1}^{\prime}&q_{2}^{\prime}\\\
p_{1}^{\prime}&p_{2}^{\prime}\end{pmatrix}\det\begin{pmatrix}q_{1}&q_{2}\\\
p_{1}&p_{2}\end{pmatrix}^{-1}=(q_{1}^{\prime}p_{2}^{\prime}-p_{1}^{\prime}q_{2}^{\prime})(q_{1}p_{2}-p_{1}q_{2})^{-1},$

which is equal to $1$ if and only if
$p_{1}^{\prime}q_{2}^{\prime}-q_{1}^{\prime}p_{2}^{\prime}=p_{1}q_{2}-q_{1}p_{2}$.

∎

We are now ready to evaluate $\mathcal{E}_{d}^{(3)}$, proving Theorem 25(iii).
Observe that the operators (Eq. (240) - (244)) appearing in our final
expression for $\mathcal{E}_{d}^{(3)}$ corresponds to the non-trivial orbits
we found in Lemma 29. For each orbit of the form
$\\{(D_{q_{1},p_{1}}D_{q_{2},p_{2}})\\}$, we have an operator that is an equal
superposition over states that are proportional to
$|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{-q_{1}-q_{2},-p_{1}-p_{2}}\rrangle$
(and recall that states of this form are the states not ruled out by Lemma
28). In this sense, $|\Phi_{1}\rrangle$ corresponds to the orbit in Lemma
28b), $|\Phi_{2}\rrangle$ corresponds to the orbit in c), $|\Phi_{3}\rrangle$
corresponds to the orbit in d) for $k=d-1$, the $|\Psi_{k}\rrangle$’s
correspond to the orbits in d) for $k\neq d-1$, and the
$|\Upsilon_{l}\rrangle$’s correspond to the orbits in e).

###### Proof of Theorem 25(iii).

Inserting resolutions of the identity (Eq. (237)) and using Lemma 28, we have

$\displaystyle\mathcal{E}_{d}^{(3)}$
$\displaystyle=\sum_{q_{1},p_{1},q_{2},p_{2},q_{1}^{\prime},p_{1}^{\prime},q_{2}^{\prime},p_{2}^{\prime}\in\mathbb{Z}_{d}}\llangle
D_{q_{1},p_{1}}|\llangle D_{q_{2},p_{2}}|\llangle
D_{-q_{1}-q_{2},-p_{1}-p_{2}}|\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{-q_{1}^{\prime}-q_{2}^{\prime},-p_{1}^{\prime}-p_{2}^{\prime}}\rrangle$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{-q_{1}-q_{2},-p_{1}-p_{2}}\rrangle\llangle
D_{q_{1}^{\prime},p_{1}^{\prime}}|\llangle
D_{q_{2}^{\prime},p_{2}^{\prime}}|\llangle
D_{-q_{1}^{\prime}-q_{2}^{\prime},-p_{1}^{\prime}-p_{2}^{\prime}}|$
$\displaystyle=\sum_{q_{1},p_{1},q_{2},p_{2},q_{1}^{\prime},p_{1}^{\prime},q_{2}^{\prime},p_{2}^{\prime}\in\mathbb{Z}_{d}}\llangle
D_{q_{1},p_{1}}|\llangle D_{q_{2},p_{2}}|\llangle
D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}|\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}\rrangle\llangle
D_{q_{1}^{\prime},p_{1}^{\prime}}|\llangle
D_{q_{2}^{\prime},p_{2}^{\prime}}|\llangle
D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}|,$

using the fact that $D_{-q_{1}-q_{2},-p_{1}-p_{2}}\propto
D_{-q_{1},-p_{1}}D_{-q_{2},-p_{2}}=D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}$
in the second equality.

Next, note from the definition of the orbits (Eq. (251)) that there exists
$U\in\mathrm{Cl}_{d}$ such that $\llangle D_{q_{1},p_{1}}|\llangle
D_{q_{2},p_{2}}|\llangle
D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}|\mathcal{U}^{\otimes
3}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle\neq
0$ if and only if $(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$ and
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ are in
the same orbit. Hence, since since $\mathcal{E}_{d}^{(3)}$ is the average of
$\mathcal{U}^{\otimes 3}$ for $U\in\mathrm{Cl}_{d}$, $\llangle
D_{q_{1},p_{1}}|\llangle D_{q_{2},p_{2}}|\llangle
D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}|\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle\neq
0$ only if $(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$ and
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ are in
the same orbit, so we can write

$\displaystyle\mathcal{E}_{d}^{(3)}$
$\displaystyle=\sum_{\text{orbits}}\sum_{\begin{subarray}{c}(D_{q_{1},p_{1}},D_{q_{2},p_{2}}),\\\
(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})\in\text{
orbit}\end{subarray}}\llangle D_{q_{1},p_{1}}|\llangle
D_{q_{2},p_{2}}|\llangle
D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}|\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\times|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}\rrangle\llangle
D_{q_{1}^{\prime},p_{1}^{\prime}}|\llangle
D_{q_{2}^{\prime},p_{2}^{\prime}}|\llangle
D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}|.$

Now, we show that the overlap $\llangle D_{q_{1},p_{1}}|\llangle
D_{q_{2},p_{2}}|\llangle
D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}|\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle$
is the same within each orbit, i.e., for any orbit, there is a constant $c$
such that

$\llangle D_{q_{1},p_{1}}|\llangle D_{q_{2},p_{2}}|\llangle
D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}|\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle=c$

for all $(D_{q_{1},p_{1}},D_{q_{2},p_{2}})$ and
$(D_{q_{1}^{\prime},p_{1}^{\prime}},D_{q_{2}^{\prime},p_{2}^{\prime}})$ in
that orbit. To see this, note from the definition of orbits that for any other
$(D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}},D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}})$
in the orbit, there exist $U\in\mathrm{Cl}_{d}$ such that
$\mathcal{U}^{\otimes
2}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle=\alpha|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}\rrangle|D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}\rrangle$
for some phase $\alpha$. Then,

$\displaystyle\mathcal{U}^{\otimes
3}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle$
$\displaystyle=\alpha|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}\rrangle|D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}\rrangle\alpha^{*}|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}^{\dagger}D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}^{\dagger}\rrangle$
$\displaystyle=|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}\rrangle|D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}\rrangle|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}^{\dagger}D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}^{\dagger}\rrangle,$
(253)

so we can use Fact 4 to write

$\displaystyle\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle$
$\displaystyle=\mathcal{E}_{d}^{(3)}\mathcal{U}^{\otimes
3}|D_{q_{1}^{\prime},p_{1}^{\prime}}\rrangle|D_{q_{2}^{\prime},p_{2}^{\prime}}\rrangle|D_{q_{1}^{\prime},p_{1}^{\prime}}^{\dagger}D_{q_{2}^{\prime},p_{2}^{\prime}}^{\dagger}\rrangle$
$\displaystyle=\mathcal{E}_{d}^{(3)}|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}\rrangle|D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}\rrangle|D_{q_{1}^{\prime\prime},p_{1}^{\prime\prime}}^{\dagger}D_{q_{2}^{\prime\prime},p_{2}^{\prime\prime}}^{\dagger}\rrangle.$

Thus, noting that the states defined in Eqs. (240) - (244) are precisely the
equal superpositions over the
$|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}\rrangle$
for each orbit, we have

$\displaystyle\mathcal{E}_{d}^{(3)}=a|I\rrangle\llangle I|^{\otimes
3}+b_{1}|\Phi_{1}\rrangle\llangle\Phi_{1}|+b_{2}|\Phi_{2}\rrangle\llangle\Phi_{2}|+b_{3}|\Phi_{3}\rrangle\llangle\Phi_{3}|+\sum\limits_{k=1}^{d-2}c_{k}|\Psi_{k}\rrangle\llangle\Psi_{k}|+\sum\limits_{l=1}^{d-1}d_{l}|\Upsilon_{l}\rrangle\llangle\Upsilon_{l}|$
(254)

for constants $a,b_{1},b_{2},b_{3},c_{k},d_{l}$. Since $\mathcal{E}_{d}^{(3)}$
is a projector (Fact 3), each of these constants must be either $0$ or $1$.
Finally, we see from Eq. (253) and the definition of orbits (Eq. (251)) that
for any $U\in\mathrm{Cl}_{d}$, $\mathcal{U}^{\otimes 3}$ permutes the basis
states of the form
$|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}\rrangle$
within each orbit, without incurring any phases. Since each of the states
appearing in Eq. (254) is an equal superposition over these
$|D_{q_{1},p_{1}}\rrangle|D_{q_{2},p_{2}}\rrangle|D_{q_{1},p_{1}}^{\dagger}D_{q_{2},p_{2}}^{\dagger}\rrangle$
states for each orbit, it follows that $\mathcal{U}^{\otimes 3}$ stabilises
each of these states, and hence so does $\mathcal{E}_{d}^{(3)}$. Therefore,
each of the constants in Eq. (254) must be equal to $1$. ∎

###  H.4 Formulae for generalized Clifford shadows

Having evaluating the $2$\- and $3$-fold twirl channels for $\mathrm{Cl}_{d}$
(Theorem 25), we can now characterise the classical shadows corresponding to
the uniform distribution over generalized Cliffords, in any prime dimension.

####  H.4.1 Measurement channel

First, we determine the classical shadows measurement channel $\mathcal{M}$,
by substituting our expression for the $2$-fold twirl channel from Theorem
25(ii) into Eq. (225).

######  Theorem 30 (Generalized Clifford shadows measurement channel).

Let $\mathcal{M}$ be the quantum channel defined in Eq. (224). When
$\mathcal{D}$ is the uniform distribution over $\mathrm{Cl}_{d}$ for some
prime $d$ and $\mathcal{B}$ is the computational basis
$\\{|j\rangle:j\in\mathbb{Z}_{d}\\}$, $\mathcal{M}$ is given by

$\mathcal{M}(A)=\frac{1}{d+1}(\mathrm{tr}(A)I+A),$ (255)

for an arbitrary operator $A$ (equivalently,
$\mathcal{M}=\frac{d}{d+1}|I\rrangle\llangle I|+\frac{1}{d+1}\mathcal{I}$).
Therefore, $\mathcal{M}$ is invertible, with inverse given by

$\mathcal{M}^{-1}(A)=(d+1)A-\mathrm{tr}(A)I.$ (256)

###### Proof.

Since $\mathop{{}\mathbb{E}}_{U\sim\mathcal{D}}\mathcal{U}^{\otimes
2}=\mathcal{E}_{d}^{(2)}$ when $\mathcal{D}$ is the uniform distribution over
$\mathrm{Cl}_{d}$, Eq. (225) becomes

$\mathcal{M}(A)=\mathrm{tr}_{1}\left[\sum_{j\in\mathbb{Z}_{d}}\mathcal{E}_{d}^{(2)}(|j\rangle\langle
j|^{\otimes 2})(A\otimes I)\right].$

We can simplify this by observing that $\mathcal{E}_{d}^{(2)}(|j\rangle\langle
j|^{\otimes 2})=\mathcal{E}_{d}^{(2)}(|0\rangle\langle 0|^{\otimes 2})$ for
any $j\in\mathbb{Z}_{d}$. This is because for any $j\in\mathbb{Z}_{d}$, there
exists $U\in\mathrm{Cl}_{d}$ (specifically, $U=(X^{j})^{\dagger}$ such that
$\mathcal{U}(|j\rangle\langle j|)=|0\rangle\langle 0|$. Hence, by Fact 4,
$\mathcal{E}_{d}^{(2)}(|j\rangle\langle j|^{\otimes
2})=\mathcal{E}_{d}^{(2)}\circ\mathcal{U}^{\otimes 2}(|j\rangle\langle
j|^{\otimes 2})=\mathcal{E}_{d}^{(2)}(|0\rangle\langle 0|^{\otimes 2})$. Thus,

$\mathcal{M}(A)=d\,\mathrm{tr}_{1}\left[\mathcal{E}_{d}^{(2)}(|0\rangle\langle
0|^{\otimes 2})(A\otimes I)\right].$ (257)

We now calculate $\mathcal{E}_{d}^{(2)}(|0\rangle\langle 0|^{\otimes 2})$
using Theorem 25(ii). From Eq. (9), it is easy to see that

$|0\rangle\langle
0|=\frac{1}{d}\sum_{z\in\mathbb{Z}_{d}}Z^{z}=\frac{1}{d}\sum_{z\in\mathbb{Z}_{d}}D_{0,z}.$

Defining $|\Pi_{0}\rrangle\equiv|0\rangle\langle 0|$ for convenience, this
becomes

$|\Pi_{0}\rrangle=\frac{1}{\sqrt{d}}\sum_{z\in\mathbb{Z}_{d}}|D_{0,z}\rrangle,$
(258)

and we want to find

$\displaystyle\mathcal{E}_{d}^{(2)}|\Pi_{0}\rrangle^{\otimes 2}$
$\displaystyle=(|I\rrangle\llangle I|^{\otimes
2}+|\Phi_{0}\rrangle\llangle\Phi_{0}|)|\Pi_{0}\rrangle^{\otimes 2}.$

From Eqs. (236) and (258), we obtain

$\displaystyle\llangle D_{q,p}|\Pi_{0}\rrangle=\llangle
D_{q,p}^{\dagger}|\Pi_{0}\rrangle$
$\displaystyle=\frac{1}{\sqrt{d}}\delta_{q,0},$ (259)

so we have

$\llangle I|^{\otimes 2}|\Pi_{0}\rrangle^{\otimes 2}=\llangle
D_{0,0}|\Pi_{0}\rrangle^{2}=\frac{1}{d}$

and

$\displaystyle\llangle\Phi_{0}|(|\Pi_{0}\rrangle^{\otimes 2})$
$\displaystyle=\frac{1}{\sqrt{d^{2}-1}}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\llangle D_{q,p}|\llangle
D_{q,p}^{\dagger}|(|\Pi_{0}\rrangle^{\otimes
2})=\frac{1}{\sqrt{d^{2}-1}}\frac{d-1}{d}.$

Thus,

$\displaystyle\mathcal{E}_{d}^{(2)}|\Pi_{0}\rrangle^{\otimes 2}$
$\displaystyle=\frac{1}{d}|I\rrangle|I\rrangle+\frac{1}{d^{2}-1}\frac{d-1}{d}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|D_{q,p}\rrangle|D_{q,p}^{\dagger}\rrangle$
$\displaystyle=\frac{1}{d^{2}}I\otimes
I+\frac{1}{d^{2}(d+1)}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}D_{q,p}\otimes D_{q,p}^{\dagger}$
$\displaystyle=\frac{1}{d(d+1)}I\otimes
I+\frac{1}{d^{2}(d+1)}\sum_{q,p\in\mathbb{Z}_{d}}D_{q,p}\otimes
D_{q,p}^{\dagger}.$

Substituting this into Eq. (257) gives

$\displaystyle\mathcal{M}(A)$
$\displaystyle=\frac{1}{d+1}\mathrm{tr}(A)I+\frac{1}{d(d+1)}\sum_{q,p\in\mathbb{Z}_{d}}\mathrm{tr}(D_{q,p}A)D_{q,p}^{\dagger}$
$\displaystyle=\frac{1}{d+1}\mathrm{tr}(A)I+\frac{1}{d+1}A,$

where we use the fact that
$\\{\frac{1}{\sqrt{d}}D_{q,p}^{\dagger}\\}_{q,p\in\mathbb{Z}_{d}}$ forms a
Hilbert-Schmidt orthonormal basis to obtain the last equality. ∎

Thus, for the uniform distribution over generalized Cliffords (in prime
dimension), $\mathcal{M}$ and $\mathcal{M}^{-1}$ are depolarising channels.
This is reminiscent of the fact that the measurement channel for the classical
shadows associated with $n$-qubit Clifford circuits, analysed in Ref. [7], is
also a depolarising channel. In fact, the depolarising parameter in that case
has the same dependence on the dimension, i.e., the measurement channel has
the same form as that in Eq. (255), but with $d$ replaced by $2^{n}$. Note,
however, that $n$-qubit Clifford circuits are different from generalized
single-qudit Cliffords with $d=2^{n}$ (except for $n=1$). Actually, it can be
shown that the measurement channel for generalized Cliffords with $d=2^{n}$ is
not a depolarising channel, though this calculation is outside the scope of
this paper.

Substituting Eq. (256) into Eq. (226), the classical shadow estimator when $D$
is the uniform distribution over $\mathrm{Cl}_{d}$ has the form

$\hat{\rho}=(d+1)\hat{U}^{\dagger}|\hat{b}\rangle\langle\hat{b}|\hat{U}-I,$

and the estimator for $\mathrm{tr}(O\rho)$ for any observable $O$ is given by

$\displaystyle\hat{o}=\mathrm{tr}(O\mathcal{M}^{-1}(\hat{U}^{\dagger}|\hat{b}\rangle\langle\hat{b}|\hat{U}))$
$\displaystyle=(d+1)\mathrm{tr}(O\hat{U}^{\dagger}|\hat{b}\rangle\langle\hat{b}|\hat{U})-\mathrm{tr}(O)$
(260)
$\displaystyle=(d+1)\langle\hat{b}|\hat{U}O\hat{U}^{\dagger}|b\rangle-\mathrm{tr}(O).$
(261)

####  H.4.2 Variance

With both $\mathcal{E}_{d}^{(3)}$ and $\mathcal{M}^{-1}$ in hand, we can now
evaluate the variance, Eq. (231), for any operator $O$.

######  Theorem 31 (Generalized Clifford shadows variance).

Let $\hat{o}$ denote the classical shadows estimator for the expectation value
$\mathrm{tr}(O\rho)$, defined as in Eq. (227). When $\mathcal{D}$ is the
uniform distribution over $\mathrm{Cl}_{d}$ for some prime $d$ and
$\mathcal{B}$ is the computational basis $\\{|j\rangle:j\in\mathbb{Z}_{d}\\}$,
the variance $\mathrm{Var}[\hat{o}]$ of $\hat{o}$ for any observable $O$ is
given by

$\mathrm{Var}[\hat{o}]=\frac{d+1}{d}\|O_{0}\|_{\text{H-S}}^{2}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}(\widetilde{D}_{q,p}^{\dagger}O_{0})\mathrm{tr}(\widetilde{D}_{kq,kp}^{\dagger}O_{0}^{\dagger})\mathrm{tr}(\widetilde{D}_{kq,kp}\widetilde{D}_{q,p}\rho)-|\mathrm{tr}(O_{0}\rho)|^{2},$
(262)

where $O_{0}\coloneqq O-\mathrm{tr}(O)I/d$ denotes the traceless part of $O$,
and $\widetilde{D}_{q,p}\coloneqq UD_{q,p}U^{\dagger}$ for any choice of
$U\in\mathrm{Cl}_{d}$.

###### Proof.

From Eq. (232), we have

$\displaystyle\mathrm{Var}[\hat{o}]$
$\displaystyle=\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]-|\mathop{{}\mathbb{E}}[\hat{o}]|^{2}$
$\displaystyle=\mathrm{tr}\left[\sum_{j\in\mathbb{Z}_{d}}\mathcal{E}_{d}^{(3)}(|j\rangle\langle
j|^{\otimes
3})\left(\mathcal{M}^{-1}(O_{0})\otimes\mathcal{M}^{-1}(O_{0}^{\dagger})\otimes\rho\right)\right]-|\mathrm{tr}(O_{0}\rho)|^{2},$

where $\mathcal{M}^{-1}$ in this case is given by Eq. (256) in Theorem 30. For
essentially the same reason that $\mathcal{E}_{d}^{(2)}(|j\rangle\langle
j|^{\otimes 2})=\mathcal{E}_{d}^{(2)}(|0\rangle\langle 0|^{\otimes 2})$ for
any $j\in\mathbb{Z}_{d}$ (see the proof of Theorem 30), we also have
$\mathcal{E}_{d}^{(3)}(|j\rangle\langle j|^{\otimes
3})=\mathcal{E}_{d}^{(3)}(|0\rangle\langle 0|^{\otimes 3})$ for any
$j\in\mathbb{Z}_{d}$. Thus, the first term becomes

$\displaystyle\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]$
$\displaystyle=d\,\mathrm{tr}\left[\mathcal{E}_{d}^{(3)}(|0\rangle\langle
0|^{\otimes
3})\left(\mathcal{M}^{-1}(O_{0})\otimes\mathcal{M}^{-1}(O_{0}^{\dagger})\otimes\rho\right)\right].$

Now, we notate $|\Pi_{0}\rrangle\equiv|0\rangle\langle 0|$,
$|\rho\rrangle\equiv\rho$, $|O\rrangle\equiv O$, $|O\rrangle^{\dagger}\equiv
O^{\dagger}$ in Liouville notation (note that here, we deviate from the
convention in Eq. (235), in that $|\rho\rrangle$, $|O\rrangle$, and
$|O\rrangle^{\dagger}$ are not necessarily normalised). Then, since
$\mathcal{E}_{d}^{(3)}(|0\rangle\langle 0|^{\otimes 3})$ is Hermitian, we can
rewrite

$\displaystyle\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]$
$\displaystyle=d\,\llangle\Pi_{0}|^{\otimes
3}\mathcal{E}_{d}^{(3)\dagger}(\mathcal{M}^{-1}|O_{0}\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle).$

We begin by calculating $\mathcal{E}_{d}^{(3)}$. By Theorem 25(iii),

$\displaystyle\mathcal{E}_{d}^{(3)}|\Pi_{0}\rrangle^{\otimes
3}=\left(|I\rrangle\llangle I|^{\otimes
3}+|\Phi_{1}\rrangle\llangle\Phi_{1}|+|\Phi_{2}\rrangle\llangle\Phi_{2}|+|\Phi_{3}\rrangle\llangle\Phi_{3}|+\sum\limits_{k=1}^{d-2}|\Psi_{k}\rrangle\llangle\Psi_{k}|+\sum\limits_{l=1}^{d-1}|\Upsilon_{l}\rrangle\llangle\Upsilon_{l}|\right)|\Pi_{0}\rangle^{\otimes
3}.$

Using Eqs. (258) and (259), we find

$\displaystyle\llangle I|^{\otimes 3}|\Pi_{0}\rrangle^{\otimes
3}=\frac{1}{d^{3/2}},$
$\displaystyle\llangle\Phi_{i}|(|\Pi_{0}\rrangle^{\otimes
3})=\frac{1}{d^{3/2}}\frac{d-1}{\sqrt{d^{2}-1}}\qquad\text{for
$i\in\\{1,2,3\\}$},$ $\displaystyle\llangle\Psi_{k}|(|\Pi_{0}\rrangle^{\otimes
3})=\frac{1}{d^{3/2}}\frac{d-1}{\sqrt{d^{2}-1}}\qquad\text{for
$k\in\\{1,\dots,d-2\\}$},$
$\displaystyle\llangle\Upsilon_{l}|(|\Pi_{0}\rrangle^{\otimes
3})=0\qquad\text{for $l\in\\{1,\dots,d-1\\}$},$

so

$\displaystyle\mathcal{E}_{d}^{(3)}|\Pi_{0}\rrangle^{\otimes
3}=\frac{1}{d^{3/2}}|I\rrangle^{\otimes
3}+\frac{1}{d^{3/2}}\frac{d-1}{\sqrt{d^{2}-1}}\left(|\Phi_{1}\rrangle+|\Phi_{2}\rrangle+|\Phi_{3}\rrangle+\sum_{k=1}^{d-2}|\Psi_{k}\rrangle\right).$

Hence,

$\displaystyle\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]=\left[\frac{1}{\sqrt{d}}\llangle
I|^{\otimes
3}+\frac{1}{\sqrt{d}}\frac{d-1}{\sqrt{d^{2}-1}}\left(\llangle\Phi_{1}|+\llangle\Phi_{2}|+\llangle\Phi_{3}|+\sum_{k=1}^{d-2}\llangle\Psi_{k}|\right)\right](\mathcal{M}^{-1}|O_{0}\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle).$

From Eq. (256), we have

$\mathcal{M}^{-1}|D_{q,p}\rrangle=(d+1)|D_{q,p}\rrangle-\delta_{q,0}\delta_{p,0}|I\rrangle.$

Using this along with the fact that $\mathcal{M}^{-1}$ is self-adjoint, i.e.,
$\llangle A|\mathcal{M}^{-1}|B\rrangle=\llangle\mathcal{M}^{-1}(A)|B\rrangle$,
we obtain

$\displaystyle\frac{1}{\sqrt{d}}\llangle I|^{\otimes
3}(\mathcal{M}^{-1}|O_{0}\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle)=\frac{1}{d^{2}}|\mathrm{tr}(O_{0})|^{2}=0,$
$\displaystyle\frac{1}{\sqrt{d}}\frac{d-1}{\sqrt{d^{2}-1}}\llangle\Phi_{1}|(\mathcal{M}^{-1}|O_{0}\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle)=\frac{1}{d}\mathrm{tr}(O_{0})\mathrm{tr}(\rho
O_{0}^{\dagger})-\frac{1}{d^{2}}|\mathrm{tr}(O_{0})|^{2}=0$
$\displaystyle\frac{1}{\sqrt{d}}\frac{d-1}{\sqrt{d^{2}-1}}\llangle\Phi_{2}|(\mathcal{M}^{-1}|O_{0}\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle)=\frac{1}{d}\mathrm{tr}(O^{\dagger})\mathrm{tr}(\rho
O_{0})-\frac{1}{d^{2}}|\mathrm{tr}(O_{0})|^{2}=0$
$\displaystyle\frac{1}{\sqrt{d}}\frac{d-1}{\sqrt{d^{2}-1}}\llangle\Phi_{3}|(\mathcal{M}^{-1}|O_{0}\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle)=\frac{d+1}{d}\mathrm{tr}(O_{0}^{\dagger}O_{0})-\frac{d+1}{d^{2}}|\mathrm{tr}(O_{0})|^{2}=\frac{d+1}{d}\|O_{0}\|_{\text{H-S}}^{2}$
$\displaystyle\frac{1}{\sqrt{d}}\frac{d-1}{\sqrt{d^{2}-1}}\llangle\Psi_{k}|(\mathcal{M}^{-1}|O\rrangle\otimes\mathcal{M}^{-1}|O_{0}^{\dagger}\rrangle\otimes|\rho\rrangle)=\frac{d+1}{d^{2}}\sum_{\begin{subarray}{c}p,q\in\mathbb{Z}_{d}\\\
(p,q)\neq(0,0)\end{subarray}}\mathrm{tr}(D_{q,p}^{\dagger}O_{0})\mathrm{tr}(D_{kq,kp}^{\dagger}O_{0}^{\dagger})\mathrm{tr}(D_{kq,kp}D_{q,p}\rho).$

Thus, we arrive at

$\displaystyle\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]$
$\displaystyle=\frac{d+1}{d}\|O_{0}\|_{\text{H-S}}^{2}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}(D_{q,p}^{\dagger}O_{0})\mathrm{tr}(D_{kq,kp}^{\dagger}O_{0}^{\dagger})\mathrm{tr}(D_{kq,kp}D_{q,p}\rho),$
(263)

giving

$\mathrm{Var}[\hat{o}]=\frac{d+1}{d}\|O_{0}\|_{\text{H-S}}^{2}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}({D}_{q,p}^{\dagger}O_{0})\mathrm{tr}({D}_{kq,kp}^{\dagger}O_{0}^{\dagger})\mathrm{tr}({D}_{kq,kp}{D}_{q,p}\rho)-|\mathrm{tr}(O_{0}\rho)|^{2}.$

To see that this actually holds with $D_{q,p}$ replaced by
$\widetilde{D}_{q,p}\coloneqq UD_{q,p}U^{\dagger}$ for any
$U\in\mathrm{Cl}_{d}$, note that since
$\mathcal{U}\circ\mathcal{E}_{d}^{(3)}\circ\mathcal{U}^{\dagger}$ for any
$U\in\mathrm{Cl}_{d}$ from Fact 4, we can replace $|\Phi_{i}\rrangle$,
$|\Psi_{k}\rrangle$, $|\Upsilon_{l}\rrangle$ with
$\mathcal{U}|\Phi_{i}\rrangle$, $\mathcal{U}|\Psi_{k}\rrangle$,
$\mathcal{U}|\Upsilon_{l}\rrangle$, respectively, in all of our calculations.
These have the same form as $|\Phi_{i}\rrangle$, $|\Psi_{k}\rrangle$,
$|\Upsilon_{l}\rrangle$ (Eqs. (240) - (244)), respectively, except with
$D_{q,p}$ replaced by $\widetilde{D}_{q,p}$, thus leading to Eq. (262). ∎

We can sanity-check this result by comparing to that for the uniform
distribution over $n$-qubit Clifford circuits; $2$ is prime, so our results
should coincide for $d=2$ (i.e., $n=1$). For $d=2$, Eq. (263) becomes

$\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]=\frac{3}{2}\|O_{0}\|^{2}_{\text{H-S}},$

since the outside sum in the second term is empty. On the other hand, For
$n$-qubit Clifford circuits, Eq. (S43) of Ref. [7] gives
$\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]=\frac{2^{n}+1}{2^{n}+2}[\mathrm{tr}(O_{0}^{2})+2\mathrm{tr}(\rho
O_{0}^{2})]$, which is

$\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]=\frac{3}{4}[\mathrm{tr}(O_{0}^{2})+2\mathrm{tr}(\rho
O_{0}^{2})]$ (264)

for $n=1$. Ref. [7] considers only Hermitian $O_{0}$, so the RHS is always
real, and $\mathrm{tr}(O_{0}^{2})=\|O_{0}\|^{2}_{\text{H-S}}$. Now note that
for $n=1$, $O_{0}$ has dimension $2$, so since it is traceless, we must have
$O_{0}^{2}=aI$ for some $a\in\mathbb{R}$. Thus,
$\|O_{0}\|^{2}_{\text{H-S}}=\mathrm{tr}(O_{0}^{2})=a\mathrm{tr}(I)=2a$, and
$\mathrm{tr}(\rho
O_{0}^{2})=a\mathrm{tr}(\rho)=a=\|O_{0}\|^{2}_{\text{H-S}}/2$, so we can write
Eq. (264)

$\mathop{{}\mathbb{E}}[|\hat{o}|^{2}]=\frac{3}{4}\left[\|O_{0}\|_{\text{H-s}}^{2}+2\frac{\|O_{0}\|_{\text{H-s}}^{2}}{2}\right]=\frac{3}{2}\|O_{0}\|_{\text{H-s}}^{2},$

which matches our result.

###  H.5 Some observables of interest

In the section, we apply our general results, Theorems 30 and 31 to some
natural observables to consider in the qudit setting.

####  H.5.1 Transition elements in stabiliser bases

First, we consider observables of the form $U|j\rangle\langle i|U^{\dagger}$
for $U\in\mathrm{Cl}_{d}$ and computational basis states $|i\rangle$,
$|j\rangle$. The expectation values of these with respect to an unknown state
$\rho$ are $\mathrm{tr}(U|j\rangle\langle i|U^{\dagger})=\langle
i|U^{\dagger}\rho U|j\rangle$, which are the quantities considered in Theorem
6.

###### Proof of Theorem 6.

Consider any generalized Clifford $U\in\mathrm{Cl}_{d}$, and any computational
basis states $|i\rangle$ and $|j\rangle$ with $i\neq j$. We have

$U|j\rangle\langle i|U^{\dagger}=UX^{i}|j-i\rangle\langle
0|(X^{i})^{\dagger}U^{\dagger}=UX^{i}|j^{\prime}\rangle\langle
0|(UX^{i})^{\dagger},$

defining $j^{\prime}\equiv j-i\enspace(\text{mod }d)$. Hence, since
$UX^{i}\in\mathrm{Cl}_{d}$, we can choose
$\widetilde{D}_{q,p}=(UX^{i})D_{q,p}(UX^{i})^{\dagger}$ in Theorem 31, to get
that for $O=U|j\rangle\langle i|U^{\dagger}$ (which is already traceless, for
$j\neq i$),

$\displaystyle\mathrm{Var}[\hat{o}]\Big{|}_{O=U|j\rangle\langle
i|U^{\dagger}}$ $\displaystyle=\frac{d+1}{d}\|U|j\rangle\langle
i|U^{\dagger}\|_{\text{H-S}}^{2}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}(\widetilde{D}_{q,p}^{\dagger}U|j\rangle\langle
i|U^{\dagger})\mathrm{tr}(\widetilde{D}_{kq,kp}^{\dagger}(U|j\rangle\langle
i|U^{\dagger})^{\dagger})\mathrm{tr}(\widetilde{D}_{kq,kp}\widetilde{D}_{q,p}\rho)$
$\displaystyle\quad-|\mathrm{tr}(U|j\rangle\langle i|U^{\dagger}\rho)|^{2}$
$\displaystyle=\frac{d+1}{d}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}(D_{q,p}^{\dagger}|j^{\prime}\rangle\langle
0|)\mathrm{tr}(D_{kq,kp}^{\dagger}|0\rangle\langle
j^{\prime}|)\mathrm{tr}(\widetilde{D}_{kq,kp}\widetilde{D}_{q,p}\rho)-|\mathrm{tr}(U|j\rangle\langle
i|U^{\dagger}\rho)|^{2}$
$\displaystyle\leq\frac{d+1}{d}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}|\langle
0|D_{q,p}^{\dagger}|j^{\prime}\rangle\langle
j^{\prime}|D_{kq,kp}^{\dagger}|0\rangle|,$

where in the second line, we use the fact that

$\mathrm{tr}(\widetilde{D}_{q,p}^{\dagger}U|j\rangle\langle
i|U^{\dagger})=\mathrm{tr}((UX^{i})D_{q,p}(UX^{i})^{\dagger}UX^{i}|j^{\prime}\rangle\langle
0|(UX_{i})^{\dagger})=\mathrm{tr}(|j^{\prime}\rangle\langle 0|)$

for our choice of $\widetilde{D}_{q,p}$, and in the third line we use
$|\mathrm{tr}(\widetilde{D}_{kq,kp}\widetilde{D}_{q,p}\rho)|\leq 1$. Now, from
Eq. (12), we have

$D_{q,p}=\omega^{qp/2}\sum_{j\in\mathbb{Z}_{d}}\omega^{jp}|j+q\rangle\langle
j|,$ (265)

from which we obtain $|\langle
0|D_{q,p}^{\dagger}|j^{\prime}\rangle=\delta_{q,j^{\prime}}$ and $|\langle
j^{\prime}|D_{kq,kp}^{\dagger}|0\rangle=\delta_{kq,-j^{\prime}}$, so

$\displaystyle\mathrm{Var}[\hat{o}]\Big{|}_{O=U|j\rangle\langle
i|U^{\dagger}}\leq\frac{d+1}{d}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\delta_{q,j^{\prime}}\delta_{kq,-j^{\prime}}=\frac{d+1}{d}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{p=1}^{d-1}\delta_{kj^{\prime},-j^{\prime}}.$

Since $i\neq j$, we have $j^{\prime}=j-i\neq 0$, so $kj^{\prime}=-j^{\prime}$
only for $k=d-1$, which is not included in the sum over $k$. Therefore, the
second term vanishes, and we have

$\mathrm{Var}[\hat{o}]\Big{|}_{O=U|j\rangle\langle
i|U^{\dagger}}\leq\frac{d+1}{d}<2,$

for any choice of $U\in\mathrm{Cl}_{d}$, and $|i\rangle,|j\rangle$ with $i\neq
j$.

Substituting this into Eq. (233), it follows that

$N_{\text{sample}}=\mathcal{O}\left(\frac{\log(M/\delta)}{\varepsilon^{2}}\right)$

single-copy measurements of $\rho$ suffice to estimate any $M$ observables of
the form $\mathrm{tr}(U|j\rangle\langle i|U^{\dagger}\rho)=\langle
i|U^{\dagger}\rho U|j\rangle$ with $U\in\mathrm{Cl}_{d}$ and $i\neq j$, within
additive error $\varepsilon$ with high probability. Since
$|\mathrm{Cl}_{d}|=\mathcal{O}(d^{5})$ [52], there are $\mathcal{O}(d^{7})$
possible such observables,
$N_{\text{sample}}=\mathcal{O}(\log(d)/\varepsilon^{2})$ suffice to estimate
all of them.

By Eq. (261), if $U^{\prime}\in\mathrm{Cl}_{d}$ was the Clifford sampled in
one iteration of the classical shadows protocol, and $|b\rangle$ was the
computational basis state that was measured, the corresponding sample of the
estimator $\hat{o}$ for $\langle i|U^{\dagger}\rho U|j\rangle$ is

$\displaystyle\mathrm{tr}(U|j\rangle\langle
i|U^{\dagger}\mathcal{M}^{-1}(U^{\prime\dagger}|b\rangle\langle
b|U^{\prime}))$ $\displaystyle=(d+1)\langle b|U^{\prime}U|j\rangle\langle
i|U^{\dagger}U^{\prime\dagger}|b\rangle.$

Since $|b\rangle,|i\rangle,|j\rangle$ are computational basis states and
$U^{\prime}U\in\mathrm{Cl}_{d}$, each of $\langle b|U^{\prime}U|j\rangle$ and
$\langle i|U^{\dagger}U^{\prime\dagger}|b\rangle$ can be efficiently computed
(classically) using the $d$-dimensional generalization of the Gottesman-Knill
theorem [53], so the classical shadows samples can be efficiently computed. ∎

On the other hand, for an observable of the form $U|0\rangle\langle
0|U^{\dagger}$ for $U\in\mathrm{Cl}_{d}$ (which subsumes all observables
$U|i\rangle\langle i|U^{\dagger}$ for $U\in\mathrm{Cl}_{d}$ and computational
states $|i\rangle$), the traceless part is $U|i\rangle\langle
i|U^{\dagger}-I/d$, and choosing $\widetilde{D}_{q,p}=UD_{q,p}U^{\dagger}$ in
Theorem 31 gives

$\displaystyle\mathrm{Var}[\hat{o}]\Big{|}_{O=U|0\rangle\langle
0|U^{\dagger}}$ $\displaystyle=\frac{d+1}{d}\left(1-\frac{1}{d}\right)$
$\displaystyle\quad+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}\left(\widetilde{D}_{q,p}^{\dagger}\Big{(}U|0\rangle\langle
0|U^{\dagger}-\frac{I}{d}\Big{)}\right)\mathrm{tr}\left(\widetilde{D}_{kq,kp}^{\dagger}\Big{(}U|0\rangle\langle
0|U^{\dagger}-\frac{I}{d}\Big{)}\right)\mathrm{tr}(\widetilde{D}_{kq,kp}\widetilde{D}_{q,p}\rho)$
$\displaystyle\quad-\left|\mathrm{tr}\left(\Big{(}U|0\rangle\langle
0|U^{\dagger}-\frac{I}{d}\Big{)}\rho\right)\right|^{2}$
$\displaystyle=1-\frac{1}{d^{2}}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\langle 0|D_{q,p}^{\dagger}|0\rangle\langle
0|D_{kq,kp}^{\dagger}|0\rangle\mathrm{tr}({D}_{kq,kp}{D}_{q,p}U^{\dagger}\rho
U)-(\langle 0|U^{\dagger}\rho U|0\rangle-1/d)^{2}.$

From Eq. (265),

$\langle 0|D_{q,p}^{\dagger}|0\rangle=\delta_{q,0},$

so the double sum becomes

$\displaystyle\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\delta_{q,0}\delta_{kq,0}\mathrm{tr}({D}_{kq,kp}{D}_{q,p}U^{\dagger}\rho
U)$
$\displaystyle=\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{p=1}^{d-1}\mathrm{tr}(D_{0,kp}D_{0,p}U^{\dagger}\rho
U)$
$\displaystyle=\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{p=1}^{d-1}\mathrm{tr}(Z^{(k+1)p}U^{\dagger}\rho
U)$
$\displaystyle=\frac{(d+1)(d-2)}{d^{2}}\mathrm{tr}\left(\sum_{p=1}^{d-1}Z^{p}U^{\dagger}\rho
U\right)$
$\displaystyle=\frac{(d+1)(d-2)}{d^{2}}\mathrm{tr}\left((d|0\rangle\langle
0|-I)U^{\dagger}\rho U\right).$

Hence,

$\displaystyle\mathrm{Var}[\hat{o}]\Big{|}_{O=U|0\rangle\langle
0|U^{\dagger}}$
$\displaystyle=1-\frac{1}{d^{2}}+\frac{(d+1)(d-2)}{d^{2}}(d\langle
0|U^{\dagger}\rho U|0\rangle-1)-(\langle 0|U^{\dagger}\rho
U|0\rangle-1/d)^{2}$ $\displaystyle=\frac{1}{d}+(d-1)\langle 0|U^{\dagger}\rho
U|0\rangle-\langle 0|U^{\dagger}\rho U|0\rangle^{2},$

which be $\Omega(d)$ for certain states $\rho$.

####  H.5.2 Displacement operators

Now, we calculate the variance for estimating the expectation values of
displacement operators $D_{a,b}$. For $(a,b)\neq(0,0)$, i.e., $D_{a,b}\neq I$,
$D_{a,b}$ is traceless, so choosing $\widetilde{D}_{q,p}=D_{q,p}$ in Theorem
31 gives

$\displaystyle\mathrm{Var}[\hat{o}]\Big{|}_{O=D_{a,b}}$
$\displaystyle=\frac{d+1}{d}\|D_{a,b}\|_{\text{H-S}}^{2}+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\mathrm{tr}(D_{q,p}^{\dagger}D_{a,b})\mathrm{tr}(D_{kq,kp}^{\dagger}D_{a,b}^{\dagger})\mathrm{tr}(D_{kq,kp}D_{q,p}\rho)-|\mathrm{tr}(D_{a,b}\rho)|^{2}$
$\displaystyle=\frac{d+1}{d}(d)+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\sum_{\begin{subarray}{c}q,p\in\mathbb{Z}_{d}\\\
(q,p)\neq(0,0)\end{subarray}}\delta_{q,a}\delta_{p,b}\delta_{kq,-a}\delta_{kp,-b}\mathrm{tr}(D_{kq,kp}D_{q,p}\rho)-|\mathrm{tr}(D_{a,b}\rho)|^{2}$
$\displaystyle=d+1+\frac{d+1}{d^{2}}\sum_{k=1}^{d-2}\delta_{ka,-a}\delta_{kb,-b}\mathrm{tr}(D_{ka,ka}D_{a,b}\rho)-|\mathrm{tr}(D_{a,b}\rho)|^{2}.$

If $a\neq 0$, then $ka=-a$ only for $k=d-1$, which does not appear in the sum
over $k$. Otherwise, we must have $b\neq 0$, and similarly no value of $k$ in
the sum satisfies $kb=-b$. Hence, for any $D_{a,b}\neq I$,

$\mathrm{Var}[\hat{o}]\Big{|}_{O=D_{a,b}}=d+1-|\mathrm{tr}(D_{a,b}\rho)|^{2}\geq
d.$

This implies that if we were to use the generalized Clifford classical shadows
to estimate $\mathrm{tr}(D_{a,b}\rho)$, and we chose the number of samples
according to the variance, via Eq. (233), the number of copies of $\rho$ would
scale with $d$ and $\varepsilon$ as $d/\varepsilon^{2}$. This is consistent
with Theorem 3. Note, however, that Theorem 3 concerns measuring
$|\mathrm{tr}(D_{a,b}\rho)|$ for all $a,b$, whereas Eq. (233) calls for
$\mathcal{O}(d/\varepsilon^{2})$ measurements even for estimating the
expectation value of just one displacement operator $D_{a,b}$.


# Crystal-like thermal transport in amorphous carbon

Jaeyun Moon  To whom correspondence should be addressed; E-mail:
<EMAIL_ADDRESS>Sibley School of Mechanical and Aerospace Engineering,
Cornell University  Zhiting Tian  Sibley School of Mechanical and Aerospace
Engineering, Cornell University

(August 28, 2024)

###### Abstract

Thermal transport properties of amorphous carbon has attracted increasing
attention due to its extreme thermal properties: It has been reported to have
among the highest thermal conductivity for bulk amorphous solids up to $\sim$
37 Wm-1K-1, comparable to crystalline sapphire ($\alpha$-Al2O3). Further,
large density dependence in thermal conductivity demonstrates a potential for
largely tunable thermal conductivity. However, mechanism behind the high
thermal conductivity and its large density dependence remains elusive due to
many variables at play. In this work, we perform large-scale ($\sim$ 105
atoms) molecular dynamics simulations utilizing a machine learning potential
based on neural networks. Through spectral decomposition of thermal
conductivity which enables a quantum correction to classical heat capacity, we
find that propagating vibrational excitations govern thermal transport in
amorphous carbon ($\sim$ 100 % of thermal conductivity) in sharp contrast to
the conventional wisdom that diffusive vibrational excitations dominate
thermal transport in amorphous solids. Instead, this remarkable behavior
resembles thermal transport in simple crystals. Moreover, our temperature
dependent spectral diffusivity and velocity current correlation analyses
reveal that the density dependent thermal conductivity originates from
anharmonicity sensitive propagating excitations. Our work suggests a novel
insight and design principle into developing mechanically hard, thermally
conductive amorphous solids.

##  I Introduction

Vibrational properties of amorphous solids are of fundamental interest due to
their anomalies compared to those of crystalline solids, including a Boson
peak in the vibrational density of states [shintani_universal_2008,
chumakov_role_2014] and an excess heat capacity at cryogenic temperatures
[zeller_thermal_1971]. Combined with small average mean free paths on the
order of interatomic distances [kittel_interpretation_1949], these anomalies
suggest different characteristics of vibrational heat carriers from those of
phonons in simple crystals. Prior seminal works have proposed to categorize
these heat carriers in amorphous solids as propagating, diffusive, and
localized vibrations depending on their transport mechanisms
[allen_thermal_1989, allen_diffusons_1999, allen_character_2003]. Propagating
vibrational excitations are similar to phonons in that they have large mean
free paths, thus transporting heat efficiently. Diffusive vibrations have ill-
defined mean free paths below $\sim$ interatomic distances and transport heat
in a random walk manner. Localized vibrations are spatially localized and are
often considered to contribute negligibly to thermal transport. These
categorizations have been used to understand various thermal properties in
numerous amorphous solids microscopically [lv_direct_2016, lv_examining_2016,
lv_non-negligible_2016, larkin_thermal_2014, moon_propagating_2018,
deangelis_thermal_2018, sorensen_heat_2020, aryana_tuning_2021]. Consensus is
that diffusive vibrations govern thermal transport in most amorphous solids
due to large degrees of disorder hindering propagation of vibrations, leading
to low thermal conductivity in dielectric amorphous solids ($\lesssim 1$
Wm-1K-1).

Recently, amorphous carbon has drawn immense attention due to its remarkable
mechanical and thermal properties [zeng_pressure_2018, lin_amorphous_2011,
shang_enhancement_2023, zhang_discovery_2022, zhang_narrow-gap_2021,
shang_ultrahard_2021]: Amorphous carbon with similar density to crystalline
diamond (3.5 g cm-3 at 300 K) can achieve larger Young’s modulus and hardness
than those of crystalline diamond [shang_ultrahard_2021] and it has among the
highest thermal conductivity reported for pure bulk amorphous solids up to
$\sim$ 37 Wm-1K-1 at 300 K, comparable to that of many crystalline solids such
as sapphire [shang_enhancement_2023, cahill_thermal_1998]. Numerous
experimental [shamsa_thermal_2006, bullen_thermal_2000, shang_elastic_2020,
shang_enhancement_2023] and computational works [lv_phonon_2016,
giri_atomic_2022, suarez-martinez_effect_2011, minamitani_relationship_2022]
have been conducted to elucidate the nature of atomic vibrations in amorphous
carbon and to explore extreme thermal properties achievable in amorphous
solids.

Prior measurements have shown that thermal conductivity has a strong positive
correlation with mass density and sp3 content, demonstrating that local atomic
environment could be an important factor in thermal transport in amorphous
carbon [shamsa_thermal_2006, bullen_thermal_2000, shang_elastic_2020,
shang_enhancement_2023]. However, thermal conductivity measurements of
amorphous carbon vary significantly in literature (from $\sim 1$ to $\sim 37$
Wm-1K-1) depending on many variables such as synthesis techniques and
conditions that lead to different densities, sample qualities, and possible
crystallinities. Accurate and systematic atomic simulation studies are, thus,
necessary to microscopically understand the thermal transport mechanisms in
amorphous carbon and these will, in larger effect, play an essential role in
understanding extreme thermal properties of amorphous solids.

Recent molecular dynamics simulations on amorphous carbon support the local
structure (density and sp3 content)-property (thermal conductivity)
relationship [lv_phonon_2016, giri_atomic_2022, suarez-martinez_effect_2011,
minamitani_relationship_2022]. By categorizing vibrations into propagating and
diffusive modes, these simulation works have further reported that diffusive
vibrations play a substantial role in thermal conduction in amorphous carbon
(up to 70 % of thermal conductivity) at 300 K, a similar behavior observed in
most other amorphous solids [moon_sub-amorphous_2016, deangelis_thermal_2018,
zhou_contribution_2017]. However, these simulations report: 1. unphysically
large vibrational density of states at high frequencies above $\sim$ 70 THz
similar to hydrogen vibrational frequencies using simple empirical potentials
[lv_phonon_2016], 2. small simulations cells, comprised of hundreds to
thousands of atoms, not able to adequately include low frequency vibrations
[lv_phonon_2016, giri_atomic_2022, minamitani_relationship_2022] 3\. classical
thermal conductivity values at room temperatures where phonon occupation is
inaccurately described below Debye temperatures ($\sim$ 2300 K for amorphous
carbon) [giri_atomic_2022, suarez-martinez_effect_2011]. Despite these
experimental and computational efforts, accurate understanding of thermal
transport in amorphous carbon is still missing.

In this work, we perform molecular dynamics simulations of amorphous carbon
using a machine learning potential based on neural networks. Large structures
made of 110,592 atoms with varying densities from 3.0 to 3.8 g cm-3 are
considered. With homogeneous non-equilibrium molecular dynamics descriptions
of thermal conductivity based on a linear response theory, we decompose
thermal conductivity into spectral contributions, enabling a quantum
correction to the phonon occupation and heat capacity. The combination of
machine learning potential with first principles accuracies, large structures,
and quantum correction to heat capacity mitigates prior limitations mentioned
before. By comparing with three widely used physical models and directly
analyzing dispersions, we show that remarkably, propagating excitations govern
thermal transport in amorphous carbon at room temperature ($\sim$ 100 % of
thermal conductivity) similar to simple crystals but unlike most amorphous
solids. This efficient mode of transport may explain the high thermal
conductivity values observed in amorphous carbon. Further, through temperature
dependent analysis, we find that large changes in density dependent thermal
conductivity observed in amorphous carbon are due to anharmonicity sensitive
propagating vibrational excitations at low frequencies below $\sim$ 10 THz.
Our results provide novel insights and design principles into a possible
mechanism behind achieving high thermal conductivity in amorphous solids.

##  II Results and discussion

###  II.1 Amorphous carbon structures

Representative amorphous solid structures at 3.0, 3.3, and 3.5 g cm-3 are
shown in Fig. 1A. Molecular dynamics simulation details are discussed in
Methods. Number of atoms having four nearest neighbors (sp3-bonded) increase
with density. A negligible number of atoms had two or five neighbors. Cutoff
distance of 1.85 Å was used to determine atomic coordination numbers
[deringer_machine_2017, sosso_understanding_2018]. A large, spatially
homogenous increase in the number of sp3-bonded atoms with increase in density
is clearly visible. All sp3-bonded atoms were subsequently identified for all
structures and their density dependent population is plotted against prior
experimental [schwan_tetrahedral_1996, fallon_properties_1993,
ferrari_density_2000, shang_elastic_2020] and computational works
[deringer_machine_2017, pastewka_describing_2008] in Fig. 1B. While all data
depict monotonic increase in sp3 population with density, some empirical
potentials such as Tersoff [deringer_machine_2017] and Brenner
[pastewka_describing_2008] underpredict the sp3 population by more than a half
compared to experimental values. Structures based on machine-learning
potentials such as Gaussian approximation potential (GAP)
[deringer_machine_2017] and NEP used here appear to have more consistent sp3
population predictions against experiments and density functional theory (DFT)
calculations. Resulting pair distribution function, $g(r)=\frac{1}{4\pi
Nnr^{2}}\sum_{i,j}\langle\delta(r-|\bm{r}_{i}-\bm{r}_{j}|)\rangle$, is
demonstrated for our 3.5 g cm-3 structure along with a prior DFT structure
[deringer_machine_2017] at the same density in Fig. 1C. Here, $N$ is the
number of atoms, $n$ is the number density, $\bm{r}_{i}$ is the atomic
position of the $i$th atom, and the angled brackets denotes an ensemble
average. As expected for an amorphous solid, broad peaks and valleys are
observed. Consistent pair distribution functions are demonstrated between the
two curves.

Figure 1: (A) Local atomic environment of tetrahedral amorphous carbon at
different mass densities. Yellow and blue atoms denote atoms with sp3 and sp2
hybridizations, respectively. At 3.0, 3.3, and 3.5 g cm-3, 58.3, 81.0, and
90.1 % of atoms are sp3-bonded, respectively. (B) Density dependent sp3
distributions of our NEP systems compared against previously reported values
for various amorphous carbon systems from atomic simulations
[deringer_machine_2017, pastewka_describing_2008] and experiments
[schwan_tetrahedral_1996, fallon_properties_1993, ferrari_density_2000,
shang_ultrahard_2021]. (C) Pair distribution function of our system with 3.5 g
cm-3 using the NEP potential (red curve) vs. prior literature using DFT
[deringer_machine_2017] (blue curve).

###  II.2 Vibrations in amorphous carbon

Using these structures, we next examine their vibrational densities of states
(DOS) using spectral velocity autocorrelations (see Fig. 2) [moon_atomic_2023,
grest_density_1981, moon_heat_2024]. Due to a large sample quantity typically
required for density of states measurements, experimental measurement
comparison is absent here. With increase in density, we observe a flatter DOS
at low frequencies below 10 THz, which signifies higher Debye sound velocities
($v_{D}$) as DOS scales as $v_{D}^{-3}$. By fitting the Debye model of
vibrational density of states below 5 THz, we obtain 12.97, 14.60, and 15.21
km s-1 for the 3.0, 3.3, and 3.5 g cm-3 structures, respectively. Vibrational
densities of states using the NEP potential appear to be more consistent with
prior DFT calculations [sosso_understanding_2018] while the DOS using a
Tersoff potential [lv_phonon_2016] predicts a significant number of vibrations
with frequencies higher than 50 THz. With the NEP potential, we then have
structural and vibrational properties with first-principles accuracies for our
large amorphous carbon structures.

Figure 2: Vibrational densities of states (DOS) of various amorphous carbon
systems. DOS for the 3.0, 3.3, and 3.5 g cm-3 structures studied here are
shown as solid blue, purple, and black curves, respectively. DFT calculations
(3.4 g cm-3) are shown as a brown dotted curve and DOS using Tersoff potential
(3.0 g cm-3) is depicted as a green dotted curve.

###  II.3 Thermal conductivity and transport mechanisms

Thermal conductivity values for all the structures were then determined by
homogeneous non-equilibrium molecular dynamics (HNEMD) [fan_homogeneous_2019,
gabourie_spectral_2021] (See Methods). Similar to the Green-Kubo (GK)
formalism [kubo_fluctuation-dissipation_1966] in equilibrium molecular
dynamics (EMD), the HNEMD method is based on a linear response theory. In
HNEMD, independent phonon calculations are not required to obtain spectral
thermal conductivity unlike other molecular dynamics based methodologies such
as Green-Kubo modal analysis (GKMA) [lv_direct_2016] and normal mode
decomposition (NMD) [mcgaughey_predicting_2014] which limit the number of
atoms due to large computational costs of constructing and diagonalizing
dynamical matrices.

Classical and quantum corrected thermal conductivity of all systems studied
here were plotted against experimental measurements and Green-Kubo thermal
conductivity independently calculated here as shown in Fig. 3. Green-Kubo
thermal conductivities were obtained by
$k=\frac{V}{3k_{B}T}\int_{0}^{\infty}dt\langle\bm{J}(t)\cdot\bm{J}(0)\rangle$
where $\bm{J}(t)$ is the heat current. Classical thermal conductivity from
HNEMD and GK thermal conductivity are within the errorbars of each other.
There have been a lot of variations in the measured thermal conductivities for
amorphous carbon depending on the synthesis techniques, thickness, and sample
quality. Nonetheless, we observe that quantum corrected thermal conductivity
for the amorphous carbon structures discussed here is on the same order of
magnitude as measurements at 300 K [bullen_thermal_2000, shamsa_thermal_2006,
chen_thermal_2000].

Figure 3: (A) Density dependent classical (red solid circles) and quantum
mechanical (black solid squares) thermal conductivity at 300 K compared
against thermal conductivity values using Green-Kubo formalism (blue solid
circles) and available measurements (green symbols). Sample synthesis
technique and thickness are: cross (filtered arc, 47 nm)
[bullen_thermal_2000], diamond (filtered cathodic vacuum arc, 18.5 to 100 nm)
[shamsa_thermal_2006], and inverse triangle (filtered cathodic vacuum arc, 20
to 100 nm) [chen_thermal_2000] (B) Classical (solid curves) and quantum
mechanical (dotted curves) spectral thermal conductivity. Shaded orange region
highlights large density dependence in spectral thermal conductivity. (C)
Thermal conductivity accumulation function of 3.0, 3.3, and 3.5 g cm-3
structures at 300 K. (D) Spectral diffusivity of 3.0, 3.3, and 3.5 g cm-3
structures at 300 K. Color schemes in (C) and (D) are the same as (B). Dashed
line representing inverse quadratic power law is a guide to the eye. Dotted
line represents diffusivity $D_{a}$ when the mean free path is interatomic
distance. Solid black line depicts a diffusivity based on a random-walk theory
$D_{RW}$.

To see the effect of the quantum correction to the heat capacity on the total
thermal conductivity in amorphous carbon, spectral thermal conductivity and
thermal conductivity accumulation functions of 3.0, 3.3, and 3.5 g cm-3
structures at 300 K are shown in Fig. 3B and 3C, respectively. It appears that
the effect of quantum correction to the specific heat on thermal conductivity
becomes prominent above $\sim$ 7 THz. This is expected as $k_{B}T$ $\sim$ 25
meV or 6 THz at room temperature. At 300 K, the quantum correction to the heat
capacity reduces the classical thermal conductivities by as much as 30 to 40
%, demonstrating the importance of considering quantum effects in these
materials. Further, it is interesting to note that density affects spectral
thermal conductivity of amorphous carbon drastically below $\sim$ 10 THz while
there is little to no effect above 10 THz. Therefore, our work reveals that
density dependent thermal conductivity in amorphous carbon widely observed in
literature at room temperature is likely due to changes in transport
properties of these low frequency vibrations.

####  II.3.1 Characterization of popagating and diffusive vibrations

We next examine spectral thermal diffusivity values for the 3.0, 3.3, and 3.5
g cm-3 structures at 300 K obtained by
$D(\omega)=\frac{k(\omega)}{C(\omega)DOS(\omega)}$ as plotted in Fig. 3D along
with some widely used diffusivity models for amorphous materials. There
currently exist a few different methods to distinguish propagating vibrations
and diffusive vibrations: some methods require solving the dynamical matrices
and subsequently examine normal mode shapes or properties (e.g. eigenvector
periodicity [seyf_method_2016]) and the other methods rely on physical models
that describe propagating and diffusive vibrations. The first route is not
computationally feasible here. Our systems ($\sim$ 105 atoms) are more than an
order of magnitude larger than typical numbers of atoms used for dynamical
matrix calculations using simple empirical potentials in literature (103 to
104 atoms) [larkin_thermal_2014, giri_atomic_2022, moon_propagating_2018].
Therefore, we resort to physical models describing propagating and diffusive
vibrations.

The first model relies on how propagating wave lifetimes depend on frequency.
At THz frequencies, lifetimes of propagating vibrations typically follow
$\omega^{-2}$. Prior works determined the frequency at which $\omega^{-2}$
dependence in lifetimes and diffusivity disappears as the transition frequency
from propagating to diffusive vibrations in amorphous solids including
amorphous carbon [larkin_thermal_2014, zhou_contribution_2017,
lv_phonon_2016]. Using this criterion, we find that the transition frequency
can be approximated to be around $\sim$ 32 THz for these amorphous carbon
structures (See Fig. S1 for more details in Supplementary Information).

The second physical model considers the lowest diffusivity ($D_{a}$) possible
by a propagating wave where group velocity is sound velocity (here, we use the
Debye velocity of $\sim$ 14 km s-1 as an average value) and mean free path is
the interatomic distance, $a$ [kittel_interpretation_1949]. The third model
considered here considers maximum diffusivity for diffusive vibrations and is
derived from a random walk theory [agne_minimum_2018] that has a linear
dependence on frequency as $D_{RW}=\frac{1}{3}\frac{\omega}{\pi n^{2/3}}$
where n is the number density. If physically realistic, these two independent
models should ideally merge at a certain frequency for amorphous carbon. These
two criteria have been previously used in amorphous solids but also in complex
crystals to determine propagating vs. diffusive normal modes that resulted in
consistent thermal conductivity with experiments [luo_vibrational_2020,
moon_propagating_2018, kim_origin_2021, cai_diffuson-dominated_2023]. Applying
these two physical models lead to nearly identical crossover frequencies at
remarkably high $\sim$ 35 THz in amorphous carbon as shown in Fig. 3D.

Three physical models considered here consistently lead to remarkably high
crossover frequencies of $\sim$ 30 to 35 THz. For these propagating
vibrational excitations then, crisp dispersions with well-defined frequencies,
group velocities, and mean free paths should be observed
[monaco_anomalous_2009, shintani_universal_2008, moon_thermal_2019,
sette_dynamics_1998, fiorentino_hydrodynamic_2023]. Representative
longitudinal and transverse dispersions for the 3.5 g cm-3 amorphous carbon
structure are shown in Fig. 4A and B (see Methods for calculation details). At
first glance, clear dispersions are demonstrated up to $\sim$ 30 THz and
$\sim$ 20 THz for longitudinal and transverse excitations, respectively. By
carefully fitting the single damped harmonic oscillator model at each
wavevector shown in Fig. S2 in the Supplementary Information, we report mean
free paths of longitudinal and transverse excitations in Fig. 4C. We find that
propagating to diffusive crossover frequencies for longitudinal and transverse
excitations are in a similar range as the other physical models at $\sim$ 30
and 35 THz, respectively. Further, the mean free paths vary by nearly four
orders of magnitude from $\sim$ 1 THz to the crossover frequencies, similar to
spectral thermal diffusivity calculated by HNEMD as shown in Fig. 3D.

Figure 4: Wavevector and frequency resolved (A) longitudinal and (B)
transverse velocity current correlations for the 3.5 g cm-3 structure.
Estimated crossover frequencies from propagating to diffusive vibrations are
noted as $\omega_{L,c}$ and $\omega_{T,c}$ for longitudinal and transverse
directions, respectively. (C) Extracted mean free paths ($\lambda$) for both
longitudinal (red circles) and transverse (blue circles) excitations.
Horizontal black line is a guide to eye representing $\lambda$ approaching
interatomic distance ($a$) of amorphous carbon.

All evidences based on the three physical models and dispersion analysis are
pointing to the consequence that despite having intrinsically disordered
structures, propagating vibrations that exist below 30 to 35 THz are
responsible for thermal transport ($\sim$ 100 % of thermal conductivity) in
amorphous carbon at 300 K as demonstraed in Fig. 3C unlike most amorphous
solids reported in literature where diffusive vibrations dominate thermal
conduction [deangelis_thermal_2018, wingert_thermal_2016]. This behavior is
remarkably similar to thermal transport in simple dielectric crystals. Our
finding is in contrast with prior molecular dynamics predictions of amorphous
carbon where a substantial contribution from diffusive vibrations is reported
(up to 70 % of thermal conductivity) [giri_atomic_2022, lv_phonon_2016]. We
attribute this difference to a combination of three possible factors: 1.
Simple empirical potentials may lack accuracy in predicting vibrational
properties [giri_atomic_2022, lv_phonon_2016, suarez-martinez_effect_2011], 2.
Prior simulations typically comprise of hundreds or thousands of atoms, not
being able to include vibrations with long wavelengths and small frequencies
adequately [giri_atomic_2022, lv_phonon_2016, minamitani_relationship_2022],
and 3. Due to the lack of spectral decomposition of thermal conductivity, a
quantum correction to the heat capacity was not possible, leading to
overpredictions of high frequency mode contributions [giri_atomic_2022].

####  II.3.2 Temperature dependent thermal transport

So far, we have considered how thermal occupation affects thermal conductivity
and the role of propagating vibrations in amorphous carbon at 300 K. We next
explore the physical nature of scattering mechanisms of these propagating
vibrations through examining temperature dependent thermal transport at 100,
300, and 500 K as shown in Fig. 5. Only the quantum corrected thermal
conductivity values are plotted in Fig. 5A along with some prior measurements
[bullen_thermal_2000, shamsa_thermal_2006]. Similar to density dependent
thermal conductivity (see Fig. 3A), some differences between our predictions
and measurements may arise from a variety of both experimental and
computational origins including deposition techniques, sample thicknesses, and
accuracy of the interatomic potential used here. A similar order of magnitude
in thermal conductivity among the measurements and predictions is,
nonetheless, shown. In all systems, a typical glass-like thermal conductivity
temperature dependence is observed: thermal conductivity increases with
increase in temperature.

Spectral diffusivity varies by $\leavevmode\nobreak\ \sim$ 6 orders of
magnitude (see Fig. 3D). Therefore, diffusivity ratios for 100 K, 300 K, and
500 K are used as our metric of temperature effect for clear visualization.
These ratios are representatively plotted for the 3.5 g cm-3 structure in Fig.
5B. Results for other amorphous carbon systems considered in this work are
generally consistent with the 3.5 g cm-3 system and are not shown for
redundancy. Above $\sim$ 10 THz, both ratios fluctuate around unity,
demonstrating that scattering mechanisms above $\sim$ 10 THz are nearly
temperature invariant and may originate from structural disorder in amorphous
carbon. Below $\sim$ 10 THz, we observe more significant temperature effects
in the spectral diffusivity. Generally, $D_{100K}/D_{500K}$ is greater than
$D_{100K}/D_{300K}$ in this frequency range, indicating that diffusivity
monotonically decreases with increase in temperature. We attribute the
temperature dependence of thermal diffusivity to anharmonic quasi-particle
interactions in this frequency range. Similar trends are observed in
independently calculated temperature dependent dispersion linewidths shown in
Fig. S3. The trend that the ratios increase as frequency decreases is thought
to originate from vibrations with larger periods (hence, larger wavelengths)
becoming less sensitive to local disorders.

Figure 5: (A) Temperature dependent thermal conductivity for 3.0 (blue solid
circles), 3.3 (purple solid circles), and 3.5 g cm-3 (black solid circles)
structures compared against prior measurements (green solid diamond: filtered
cathodic vacuum arc, 3.3 g cm-3 [shamsa_thermal_2006] and green cross:
filtered arc, 2.8 g cm-3 [bullen_thermal_2000]). These data and density
dependent thermal conductivity data with the same symbols are from the same
prior studies. (B) Spectral diffusivity ratios, $D_{100K}/D_{500K}$ (red open
circles) and $D_{100K}/D_{300K}$ (blue open circles). Diffusivity ratio of
unity is marked as black holizontal line as a guide to the eye.

Our findings of the transition of scattering mechanisms from anharmonicity
sensitive regime to disorder dominated regime are consistent with prior
experimental works on various amorphous solids ranging from network systems
[kim_origin_2021, baldi_anharmonic_2014, masciovecchio_evidence_2006,
ruocco_high-frequency_2001, ruffle_glass-specific_2006] to polymers
[sette_collective_1995, masciovecchio_inelastic_2004]: At low frequencies
below $\sim$ 10 to 100 GHz, picosecond acoustics and Brillouin light
scattering have shown strong temperature dependent mean free paths of acoustic
excitations while negligible temperature dependence in the inelastic peak
widths above peak frequency of 1 THz from inelastic X-ray and neutron
scattering measurements have been reported. The high transition frequency of
scattering mechanisms at $\sim$ 10 THz observed here in amorphous carbon may
be due to relatively high sound velocities [moon_thermal_2019] from strong C-C
bonds. At a given wavevector, higher frequency acoustic waves are supported in
amorphous carbon vs. those with low sound velocities.

Based on our density and temperature dependent vibrational and thermal
characterizations, a novel insight into how density affects microscopic
thermal transport in amorphous carbon emerges. It was observed in Fig. 2 that
the enhancement in thermal conductivity with increase in density in amorphous
carbon is mostly due to low frequency propagating vibrations below $\sim$ 10
THz. From our temperature dependent diffusivity analysis, we further conclude
that more specifically, anharmonicity sensitive propagating vibrations are
responsible for variations in thermal conductivity in amorphous carbon due to
different densities studied here. We expect that our findings here can be
extended to other amorphous solids, especially those with local tetrahedral
orders including amorphous silicon, germanium, and silica. Direct
verifications of our predictions may be possible through inelastic X-ray or
neutron scattering measurements of vibrational density of states and
dispersion of amorphous carbon.

##  III Conclusions

Using a machine learning potential based on neural networks for large-scale
amorphous carbon structures in molecular dynamics, we elucidate the
microscopic mechanism behind high thermal conductivity values in amorphous
carbon, among the highest in amorphous solids. Through spectral decomposition
of thermal conductivity, quantum correction to the heat capacity, and
dispersion analysis, we demonstrate that propagating vibrational excitations,
present up to 30 $\sim$ 35 THz, dominate thermal transport in amorphous carbon
at room temperature ($\sim$ 100 % of thermal conductivity). Further,
temperature dependent analysis of thermal diffusivity and dispersion
linewidths indicate that the origin of density dependent thermal conductivity
is due to anharmonicity sensitive propagating vibrational excitations at low
frequencies below $\sim$ 10 THz, making connections between structure and
property. Engineering anharmonicity sensitive vibrations could be key to
achieving amorphous solids with high thermal conductivity. We have suggested
an experimental approach based on inelastic scattering to test these
predictions.

##  IV Methods

###  IV.1 Molecular dynamics simulations

Crystal diamond of 110,592 atoms at various densities from 3.0 to 3.8 g cm-3
were melted at 15,000 K for 30 ps, followed by quenching to an amorphous state
at 1000 K at 100 K ps-1 in NVT ensemble (constant number of atoms, volume, and
temperature) using Graphics Processing Units Molecular Dynamics (GPUMD)
[fan_GPUMD_2022]. Resulting structures were then annealed for 1 ns at 1000 K
and were subsequently cooled to 300 K at 100 K ps-1 in NVT. Using these seed
structures, molecular dynamics outputs were recorded for 2 ns after initial
equilibration of 1 ns at desired temperatures (100, 300, and 500 K) under NVT.
Our large structures result in more statistics of propagating vibrational
excitations compared to prior molecular dynamics simulations with thousands of
atoms [lv_phonon_2016, giri_atomic_2022]. A timestep of 1 fs was used. Three
runs with different initial velocities were carried out for better statistics
for thermal and vibrational property calculations. Neuro-evolution potential
(NEP), a neural network based machine-learning potential, was utilized to
describe interatomic interactions with first-principles accuracies
[fan_GPUMD_2022, fan_neuroevolution_2021].

###  IV.2 Thermal conductivity calculations

In HNEMD, one applies an external driving force perturbation and examines the
materials response (heat current) to it as

$\bm{F}_{i}^{ext}=\bm{F}_{e}\cdot\bm{W}_{i}$ (1)

where $\bm{F}_{e}$ is the driving force parameter with the dimension of
inverse length, $\bm{W}_{i}=\sum_{j\neq i}\frac{\partial
U_{j}}{\partial\bm{r}_{ij}}\otimes\bm{r}_{ij}$ is the 3 by 3 virial tensor of
an atom, $\otimes$ represents a tensor product and $U_{j}$ is the atomic
potential energy. A small value of $\bm{F}_{e}$ at $2\times 10^{-4}$ Å-1 was
used to ensure the system is in the linear response regime.

Instantaneous heat current is related to the virial tensor by

$\bm{J}(t)=\sum_{i}\bm{W}_{i}\cdot\bm{v}_{i}.$ (2)

Ensemble average of the heat current is then directly related to the driving
force parameter and thermal conductivity by

$\langle J^{\alpha}\rangle=TV\sum_{\beta}k^{\alpha\beta}F_{e}^{\beta}$ (3)

where $\alpha$ and $\beta$ represent Cartesian directions, $T$ is temperature,
and $V$ is system volume. Thermal conductivity can therefore be determined
from Eq. 1-3. For isotropic systems like amorphous carbon, isotropic thermal
conductivity is expected. Therefore, we generally refer to thermal
conductivity as $k$. We can further decompose thermal conductivity into
spectral contributions through virial-velocity function
$Q(t)=\sum_{i}\langle\bm{W}_{i}(0)\cdot\bm{v}_{i}(t)\rangle$ as

$k(\omega)=\frac{2}{VTF_{e}}\int_{-\infty}^{\infty}Q(t)e^{i\omega t}dt.$ (4)

Due to the classical nature of molecular dynamics simulations, atomic dynamics
follows the Maxwell-Boltzmann distribution. With spectral thermal
conductivity, one could, therefore, apply a spectral ‘quantum correction’ to
the occupation number and heat capacity such that

$k_{Q}(\omega)=k(\omega)\frac{x^{2}e^{x}}{(e^{x}-1)^{2}}.$ (5)

where $x=\frac{\hbar\omega}{k_{B}T}$, $\hbar$ is the reduced Planck constant,
and $k_{B}$ is the Boltzmann constant [wang_quantum-corrected_2023]. HNEMD and
spectral decompositions of thermal conductivity have been utilized in various
materials including graphene [gabourie_spectral_2021] and amorphous silicon
[wang_quantum-corrected_2023].

###  IV.3 Dispersion analysis

Longitudinal and transverse dispersions from wavevector and frequency resolved
velocity current correlations for amorphous carbon were obtained by
[monaco_anomalous_2009, shintani_universal_2008, moon_propagating_2018,
moon_thermal_2019]

$C_{L,T}(q,\omega)=\frac{q^{2}}{2\pi\omega^{2}N}\int
dt\big{\langle}\bm{j}_{L,T}(q,t)\cdot\bm{j}_{L,T}(-q,0)\big{\rangle}e^{i\omega
t}$ (6)

where subscripts, $L,T$, refer to longitudinal or transverse branches, $q$ is
the wavevector, $\omega$ is the radial frequency, and $\bm{j}_{L,T}(q,t)$ is
given by

$\bm{j}_{L}(q,t)=\sum_{i}^{N}\big{(}\bm{v}_{i}(t)\cdot\hat{\bm{q}}\big{)}\hat{\bm{q}}e^{i\bm{q}\cdot\bm{r}_{i}(t)}$
(7)
$\bm{j}_{T}(q,t)=\sum_{i}^{N}\big{(}\bm{v}_{i}(t)-\big{(}\bm{v}_{i}(t)\cdot\hat{\bm{q}}\big{)}\hat{\bm{q}}\big{)}e^{i\bm{q}\cdot\bm{r}_{i}(t)}$
(8)

with $\bm{v}_{i}(t)$, $\bm{r}_{i}(t)$, and $\hat{q}$ representing atomic
velocities, positions, and unit wavevector, respectively.

##  V Acknowledgement

Z.T. acknowledges the support of the Department of the Navy, Office of Naval
Research under ONR award number N00014-22-1-2357. We are grateful for
discussions with Dr. Lucas Lindsay. This work used the Extreme Science and
Engineering Discovery Environment (XSEDE) Expanse under Allocation No. TG-
MAT200012. This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory,
operated under Contract No. DE-AC02-05CH11231 using NERSC award BES-
ERCAP0023621.

##  VI Author contributions

J.M. conceived the research. J.M. performed simulations and data analysis.
J.M. and Z.T. interpreted the results. J.M. wrote the paper with contributions
from Z.T.

##  VII Competing interests

The authors declare no competing interests.