# The hardness of quantum spin dynamics

Chae-Yeun Park<EMAIL_ADDRESS>Xanadu, Toronto, ON, M5G 2C8, Canada
Pablo A M Casares  Xanadu, Toronto, ON, M5G 2C8, Canada  Juan Miguel Arrazola
Xanadu, Toronto, ON, M5G 2C8, Canada  Joonsuk Huh<EMAIL_ADDRESS>Xanadu,
Toronto, ON, M5G 2C8, Canada  Department of Chemistry, Sungkyunkwan
University, Suwon 16419, Korea  SKKU Advanced Institute of Nanotechnology
(SAINT), Sungkyunkwan University, Suwon 16419, Korea  Institute of Quantum
Biophysics, Sungkyunkwan University, Suwon 16419, Korea

###### Abstract

Recent experiments demonstrated quantum computational advantage in random
circuit sampling and Gaussian boson sampling. However, it is unclear whether
these experiments can lead to practical applications even after considerable
research effort. On the other hand, simulating the quantum coherent dynamics
of interacting spins has been considered as a potential first useful
application of quantum computers, providing a possible quantum advantage.
Despite evidence that simulating the dynamics of hundreds of interacting spins
is challenging for classical computers, concrete proof is yet to emerge. We
address this problem by proving that sampling from the output distribution
generated by a wide class of quantum spin Hamiltonians is a hard problem for
classical computers. Our proof is based on the Taylor series of the output
probability, which contains the permanent of a matrix as a coefficient when
bipartite spin interactions are considered. We devise a classical algorithm
that extracts the coefficient using an oracle estimating the output
probability. Since calculating the permanent is #P-hard, such an oracle does
not exist unless the polynomial hierarchy collapses. With an anticoncentration
conjecture, the hardness of the sampling task is also proven. Based on our
proof, we estimate that an instance involving about 200 spins will be
challenging for classical devices but feasible for intermediate-scale quantum
computers with fault-tolerant qubits.

##  I Introduction

There is great interest in demonstrating that quantum computers can outperform
powerful supercomputers. The most successful paradigm so far involves showing
that sampling from the output distribution of certain quantum circuits is
classically difficult [1]. Three prominent examples of quantum sampling
problems are BosonSampling [2], random circuit sampling (RCS) [3], and
instantaneous quantum polynomial-time circuits (IQP) [4]. Theoretical works
assert that a classical computer can only efficiently simulate the same tasks
if the polynomial hierarchy (PH) collapses to the third level. Hence, under
the widely believed conjecture that PH is infinite, those sampling tasks are
classically intractable.

Based on these theoretical suggestions, several experiments have realized RCS
and a Gaussian variant of BosonSampling using noisy qubits [5] or photonic
modes [6, 7]. However, the main drawback of these quantum advantage results is
that they have little to do with practical applications. These circuits are
designed to generate nearly random bitstrings, and despite several efforts
[8], we have yet to find convincing practical applications.

On the other hand, solving quantum many-body problems is expected to be the
first useful application of quantum computers [9, 10]. Quantum many-body
Hamiltonians describing interacting quantum particles are essential in
physics, serving as models for various ranges of systems, from nuclei to
materials [11]. As we are often interested in low-temperature physics, finding
the properties of the lowest energy state, i.e., the ground state, of quantum
many-body Hamiltonians has been the central problem over decades [12]. Indeed,
many important questions in chemistry and materials science, such as the
electronic, magnetic, and mechanical properties of molecules and materials,
are answered by attributes of the ground state of these systems [13, 14].

Due to its importance, significant effort has been devoted to computing
properties of the ground state, such as the ground-state energy [15]. However,
it is well agreed that solving the ground-state energy problem of a large
quantum many-body Hamiltonian with hundreds of interacting particles is highly
challenging to a classical computer [16]. This fact is further supported by
computational complexity theory, where it has been shown that finding the
ground-state energy is QMA-hard [17] — a quantum version of NP-hard. Similar
to NP-hard problems that are regarded as difficult for classical computers, it
is widely believed that QMA-hard problems are challenging even for quantum
computers [18]. Hence, we do not expect an efficient classical algorithm for
finding the ground-state energy of an arbitrary Hamiltonian to exist.

Recently, the quantum coherent dynamics of many-body systems have also gained
considerable attention as they are related to experiments with engineered
quantum systems [19] and also provide an avenue to understand quantum chaos
[20]. The dynamical problem can, in principle, be more complicated than the
ground-state problem because it involves a manifold containing many
eigenstates, including the ground state. Despite such evidence, the hardness
of quantum dynamics has only been proven for a limited class of Hamiltonians,
such as Ising variants with commuting terms [21, 22, 23, 24], or a spin model
that can be translated into interacting bosons [25]. The hardness proof of
quantum dynamics for a broader class of Hamiltonians would be significant
given that simulating the dynamics of quantum systems is a potential first
practical application of quantum computing hardware, both for noisy [26] and
early fault-tolerant quantum computers [10]. A hardness proof can also guide
the type of experiment that should be tackled to claim quantum computational
advantage. A recent experiment [27] and follow-up classical simulation results
[28, 29, 30] demonstrate that lots of subtleties exist in proving quantum
computational advantage and why a concrete theoretical guide is vital.

In this work, we prove the classical hardness of simulating quantum spin
dynamics. Our target Hamiltonians include the Ising, XX, and Heisenberg
models, some of the most widely studied many-particle spin-$1/2$ Hamiltonians.
While some variants of IQP circuits already imply that computing the output
distribution at a certain time is classically hard for some Ising-type
Hamiltonians [21, 24, 22, 23], these results are still far from physically
interesting scenarios as the complexity is only proved at a single time point
[21, 22, 23] or for an exponentially short time [24]. The dynamics of those
Hamiltonians are also rather simple due to their commuting structure.

On the other hand, Peropadre et al. [25] studied the classical hardness of
sampling from the output distributions of the XY model, whose terms are non-
commuting. Nevertheless, their hardness result is only valid in the limit
where spins approximate hard-core bosons. This implies that their protocol
would require hundreds of thousands of spins to claim the quantum
computational advantage. In contrast, we prove that approximating the output
distribution is classically hard for a wide class of non-commuting spin
Hamiltonians without taking such a limit. We also argue that our protocol only
requires about 200 spins to demonstrate the advantage, which is more
experimentally feasible.

We utilize recently developed techniques for bounding errors from truncated
Taylor series [31, 32, 33, 34] to prove hardness. Precisely, we show that an
efficient classical algorithm estimating a polynomial function within a small
additive error can also estimate any of its high-order derivatives. This
result is particularly interesting for a certain class of Hamiltonians where
the high-order derivative of the output distribution is given by the permanent
of a matrix. Under the well-established hardness result that computing
permanent of random Gaussian matrices is #P-hard [2], we prove that estimating
the output distribution for these Hamiltonians is also classically
intractable.

With an anticoncentration conjecture, we also prove that an efficient
classical algorithm that can sample from the output distribution does not
exist unless the PH collapses. The main technical challenge here is that,
unlike other known sampling hardness results, our problem does not have the
hiding property. Roughly speaking, the hiding property refers to the fact that
the output bitstring can be hidden by randomizing the circuit instances, i.e.,
a property of the output distribution averaged over random circuit instances
does not depend on the output bitstrings [2, 4, 3]. We accommodate this
problem by proving an average-case hardness result that works even when the
output distribution anticoncentrates only for a small fraction of outcomes, at
the expense of a smaller $L_{1}$ distance between the true distribution and
the distribution we sample from. We provide numerical evidence that our models
can show such a relaxed anticoncentration property, completing arguments for
sampling hardness.

We further discuss an experimental implementation of our protocol based on
fault-tolerant intermediate-scale quantum (ISQ) hardware [35]. By implementing
the time-evolution operator using the Trotter decomposition, we roughly
estimate that a circuit with around two hundred logical qubits and less than a
billion gates is sufficient to demonstrate quantum advantage using our
protocol, which could be feasible for the first-generation of fault-tolerant
quantum computers [36, 10, 9, 37, 38].

##  II Preliminaries

Figure 1:  A pictorial description of the Heisenberg model Hamiltonian
($H_{4}$) in Eq. (4) with bipartite networks between spin groups $\sigma$ and
$\tau$ composed of $n$ spins. $J_{ij}$ denotes a random coupling strength
between $\sigma^{(i)}$ and $\tau^{(j)}$ drawn from a normal distribution.
Coupling strengths are denoted by the thickness of the corresponding edges.

Throughout the paper, we consider spin-$1/2$ models defined over a bipartite
graph, i.e., a graph with two sets of vertices $U$ and $V$ such that there is
no edge between the vertices within each set. We further restrict $U$ and $V$
to have the same number of vertices, i.e., $|U|=|V|=n$. Thus, our system
comprises a total of $2n$ spins. We denote the spin variables, i.e., Pauli
operators, for the $i$-th ($j$-th) spin of $U$ ($V$) by $\sigma^{(i)}_{x,y,z}$
($\tau^{(j)}_{x,y,z}$), where the subscript $x,y,z$ refers to the type of
Pauli operator. We consider four types of Hamiltonians defined as follows:

$\displaystyle H_{1}$
$\displaystyle=\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\sigma^{(i)}_{x}\tau^{(j)}_{x},$
(1) $\displaystyle H_{2}$
$\displaystyle=\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\bigl{[}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{z}\tau^{(j)}_{z}\bigr{]},$
(2) $\displaystyle H_{3}$
$\displaystyle=\sum_{i,j=1}^{n}\frac{J_{ij}}{2n}\bigl{[}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{y}\tau^{(j)}_{y}\bigr{]},$
(3) $\displaystyle H_{4}$
$\displaystyle=\sum_{i,j=1}^{n}\frac{J_{ij}}{2n}\bigl{[}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{y}\tau^{(j)}_{y}+\sigma^{(i)}_{z}\tau^{(j)}_{z}\bigr{]}.$
(4)

Hamiltonians $H_{1}$, $H_{3}$, and $H_{4}$ are often called the Ising, the XX,
and the Heisenberg models, respectively. See Fig. 1 for an illustration of the
Heisenberg model Hamiltonian. We will show that computing the outcome
probability distributions from the unitary evolution generated by these
Hamiltonians is classically hard.

We additionally define the following notations for binary vectors for
convenience:

  * •

We use bold letters to denote binary vectors in $\mathbb{Z}_{2}^{2n}$. For
example, $\mathbf{x}=\\{x_{1},\cdots,x_{2n}\\}$ is a vector and
$\ket{\mathbf{x}}=\ket{x_{1}}\cdots\ket{x_{2n}}$ is a product state in the
computational basis.

  * •

We denote the first and last $n$ bits of $\mathbf{x}\in\mathbb{Z}_{2}^{2n}$ by
$\mathbf{x}^{\sigma}:=\\{x_{1},\cdots x_{n}\\}$ and
$\mathbf{x}^{\tau}:=\\{x_{n+1},\cdots x_{2n}\\}$, respectively.

  * •

The Hamming weight of a vector $\operatorname{wt}(\mathbf{x})$ is the number
of 1s in the vector.

  * •

We denote a set of binary vectors with
$\operatorname{wt}(\mathbf{x}^{\sigma})=m$ and
$\operatorname{wt}(\mathbf{x}^{\tau})=n-m$ by $X_{m}$, i.e.,
$X_{m}=\\{\mathbf{x}\in\mathbb{Z}_{2}^{2n}:\operatorname{wt}(\mathbf{x}^{\sigma})=m\text{
and }\operatorname{wt}(\mathbf{x}^{\tau})=n-m\\}$.

  * •

$X:=\\{\mathbf{x}\in\mathbb{Z}_{2}^{2n}:\operatorname{wt}(\mathbf{x})=n\\}$ is
a set of binary vectors with Hamming weight $n$. Note that
$X=\cup_{i=0}^{n}X_{i}$.

We also utilize big $O$ and related notations, as summarized in Appendix A.

Quantum Spin Dynamics
$p(\mathbf{x};t)=|\langle\mathbf{x}|e^{-\mathrm{i}Ht}|\mathbf{y}_{0}\rangle|^{2}$
Conjecture 1: Estimating $p(\mathbf{x};t)$ within an additive error of
$\varepsilon 2^{-2n}$ (or $\varepsilon{2n\choose n}^{-1}$) is #P-hard. Theorem
1: Existence of classical efficient algorithms for estimating
$p(\mathbf{x};t)$ implies the collapse of the PH (by Lemma 1, 2, and 3).
Conjecture [Aaronson et al. 2011]: Estimating $\mathrm{Per}(J)^{2}$ within an
additive error of $\varepsilon n!$ is #P-hard. Conjecture 2 and Numerical
evidence: Anticoncentration of the output distribution (Paley-Zygmund
inequality). Theorem 2: Existence of an efficient classical sampler for
$p^{\prime}(\mathbf{x})$ having a small $L_{1}$ distance from $p(\mathbf{x})$
implies estimating $p(\mathbf{x})$ within an additive error of $\varepsilon
2^{-2n}/2+1/\mathrm{poly}(n)p(\mathbf{x})$. Theorem [Stockmeyer 1985]:
Existence of an efficient classical sampler implies
$\mathsf{FBPP}^{\mathsf{NP}}$ algorithm for estimating $p(\mathbf{x})$ in a
multiplicative error. Corollary 1: Existence of a classical sampler for
$p(\mathbf{x};t)$ within a polynomial small additive error collapses the PH.
Figure 2: Diagrammatic overview of relations between conjectures and theorems.
Boxes filled in gray indicate statements not developed in this study.

##  III Theoretical results and techniques

In this section, we introduce the classical hardness of approximately
estimating the output probability of quantum spin dynamics generated by four
types of Hamiltonians, $H_{1,2,3,4}$, (Theorem 1) and of sampling from the
corresponding output probability distribution (Corollary 1).

The central observation for our hardness proof is that the leading coefficient
of the Taylor expansion of the output probability contains the permanent of a
matrix. By utilizing recently developed robust polynomial regression
techniques [31, 32, 33, 34], we show that the efficient average-case
estimation of the output probabilities implies the estimation of the
permanent, which is #P-hard. In addition, with the aid of an anticoncentration
conjecture, we argue that the corresponding sampling task is also classically
challenging. The main conjectures and theorems (see Fig. 2) are followed by
evidence and proofs using the listed mathematical tools.

###  III.1 Classical hardness of computing the output distribution

We introduce a key conjecture below, which states that estimating the output
probability of quantum spin dynamics within an exponentially small additive
error is #P-hard. Many previous studies have conjectured that approximating
the output probability is #P-hard on average for classes of quantum circuits
[2, 4, 3, 22]. Our conjecture extends these arguments to quantum spin
dynamics. This conjecture is supported by one of our main results, Theorem 1,
presented later in this section.

######  Conjecture 1.

Let
$p_{1,2,3,4}(\mathbf{x};J;t):=|\braket{\mathbf{x}}{e^{-\mathrm{i}H_{1,2,3,4}t}}{\mathbf{y}_{0}}|^{2}$
where $\ket{\mathbf{y}_{0}}=\ket{\mathbf{0}}^{\sigma}\ket{\mathbf{1}}^{\tau}$,
i.e., we can write
$\mathbf{y}_{0}=\mathbf{0}^{\sigma}\mathbf{1}^{\tau}=0\cdots 01\cdots 1$. For
$\mathbf{x}\in X_{m}$ with $m\geq\sqrt{n}$ and $t_{0}=\Theta(\log n)$, the
following statements hold.

  * •

_Problem I :_ It is _#P-hard_ to approximate $p_{1}(\mathbf{x};J;t_{0})$ or
$p_{2}(\mathbf{x};J;t_{0})$ within an additive error of $\varepsilon 2^{-2n}$
with any constant probability of $p>0$ over $J\sim\mathcal{N}(0,1)^{n\times
n}$.

  * •

_Problem II :_ It is _#P-hard_ to approximate $p_{3}(\mathbf{x};J;t_{0})$ or
$p_{4}(\mathbf{x};J;t_{0})$ within an additive error of $\varepsilon{2n\choose
n}^{-1}$ with any constant probability of $p>0$ over
$J\sim\mathcal{N}(0,1)^{n\times n}$.

We provide the following three remarks in order. First, the Hamming weight $m$
must be large for $\mathbf{x}\in X_{m}$. An intuitive way to understand this
condition is that a large value of $m$ ensures that $p(\mathbf{x};J;t)$ cannot
be efficiently computed using a low-degree expansion at $t=0$. This is because
the overlap, $\braket{\mathbf{x}}{H^{k}}{\mathbf{y}_{0}}$, is non-zero only
when $k\geq m$. Thus, if one wants to compute $p(\mathbf{x};J;t)$ using the
polynomial expansion in $t$, e.g., using the Krylov subspace method [39],
$p(\mathbf{x};J;t)$ must be expanded to an order that is sufficiently larger
than $m$; this is classically difficult.

Second, the difference between the required errors $\varepsilon 2^{-2n}$ and
$\varepsilon{2n\choose n}^{-1}$ originates from the fact that $H_{3}$ and
$H_{4}$ have a global $U(1)$ symmetry, but $H_{1}$ and $H_{2}$ do not. Given
that there exist $2n$ spins, and the initial state resides in the eigenspace
of $J_{z}:=\sum_{i}\sigma^{(i)}_{z}+\sum_{j}\tau^{(j)}_{z}$ with a
corresponding eigenvalue of zero, we know that the final state is also within
this subspace since $[J_{z},H_{3,4}]=0$. As the dimension of the subspace is
$D_{II}:=|X|={2n\choose n}$, we have $\mathbb{E}_{\mathbf{x}\in
X}[p(\mathbf{x};J;t)]={2n\choose n}^{-1}$. In contrast, the dimension $D_{I}$
of the allowed space is $2^{2n}$ for $H_{1}$ and $H_{2}$, which do not have
$U(1)$ symmetry. In short, the conjecture states that approximating
$p(\mathbf{x};J;t)$ within a small factor $\varepsilon$, on an average over
$\mathbf{x}$, is classically hard.

Lastly, we require $t=\Theta(\log n)$ as the output distribution
anticoncentrates for such $t$ (throughout the paper, ‘$\log$’ stands for the
natural logarithm, unless otherwise noted). The anticoncentration property is
characterized by the fact that the output distribution $p(\mathbf{x};J;t)$
only fluctuates near (up to a constant factor) its averaged value, i.e.,
$2^{-2n}$ and ${2n\choose n}^{-1}$ for $H_{1,2}$ and $H_{3,4}$, respectively.
While the anticoncentration by itself is not a sufficient condition for
classical hardness, it supports the average-case hardness in the sense that a
trivial estimation of $p(\mathbf{x};J;t)=0$ must be incorrect [1]. Numerical
evidence for anticoncentration is provided in Sec. IV.

We now present our main theorem supporting Conjecture 1. Most importantly, we
show the classical hardness results given that the estimation of permanents of
Gaussian matrices is #P-hard as in BosonSampling [2].

######  Theorem 1.

Suppose that $J\sim\mathcal{N}(0,1)^{n\times n}$ is a randomly sampled
Gaussian matrix, $H$ is one of $H_{1,2,3,4}$, and
$p(\mathbf{x};J;t):=|\braket{\mathbf{x}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$.
Then for a given $\mathbf{x}\in X_{m}$ with $m\geq\sqrt{n}$, $t_{0}=O(\log
n)$, and a constant $0<\Delta<1$, there exists a positive constant $c_{1}$
such that the following statement holds: If there exists an efficient
classical algorithm that estimates $p(\mathbf{x};J;t_{0})$ within an additive
error of $n^{-c_{1}n^{3}\log n}$ with a probability of $\eta\geq 3/4$ over
$J$, then $\mathsf{FBPP}=\mathsf{\\#P}$.

Theorem 1 proves a weaker version of Conjecture 1. An oracle formulation of
Conjecture 1 shows the difference between the two. For example, Problem I of
Conjecture 1 is equivalent to the following statement: “If there exists an
efficient classical algorithm that estimates $p(\mathbf{x};J;t_{0})$ within an
additive error of $\varepsilon 2^{-2n}$, FP=#P.” Compared to Conjecture 1, the
theorem demands a stronger oracle for $p(\mathbf{x};J;t_{0})$, in the sense
that it requires a smaller error $n^{-c_{1}n^{3}\log n}$ and a larger
probability ($\geq 3/4$) over $J$. It also requires a slightly stronger
relation between the complexity classes for the hardness that
$\mathsf{FBPP}\neq\mathsf{\\#P}$, which is satisfied when the PH does not
collapse to the second level [40] (still widely believed to be true).

For the proof of Theorem 1, we utilize a polynomial expansion of
$p(\mathbf{x};t)=|\braket{\mathbf{x}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$
for $\mathbf{x}\in X_{m}$ (where we exclude $J$ from the arguments when they
are obvious). The expansion contains the permanent of $\tilde{J}$, a $m\times
m$ submatrix of $J$. The permanent of a matrix $A=(a_{ij})$ is defined as
$\mathrm{Per}(A):=\sum_{s\in S_{n}}\prod_{i=1}^{n}a_{i,s(i)}$ where $S_{n}$ is
the symmetric group with $n$ elements. Namely, the polynomial expansion of
$p(\mathbf{x};t)$ is given by

$\displaystyle
p(\mathbf{x};t)=\frac{\mathrm{Per}(\tilde{J})^{2}}{n^{2m}}t^{2m}+a_{2m+1}t^{2m+1}+\cdots.$
(5)

We provide a detailed derivation of this property in Appendix C. Given that
estimating the permanent of a Gaussian random matrix is #P-hard, estimating
the $2m$-th coefficient of the polynomial expansion (i.e., the coefficient for
$t^{2m}$) is a classically hard problem for $m\geq\sqrt{n}$.

We then devise an FBPP algorithm that estimates the coefficient for $t^{2m}$
when an oracle estimating $p(x;J;t)$ near $t=t_{0}$ is given. Our algorithm is
based on recently developed robust polynomial regression techniques [41, 32,
42, 31, 33, 34, 43]. We combine our algorithm with the fact that an oracle
that estimates $p(\mathbf{x};J;t_{0})$ with probability $\eta$ over
$J\sim\mathcal{N}(0,1)^{n\times n}$ can also estimate $p(\mathbf{x};J;t)$ with
probability $\eta-(3/2)n\sqrt{|(t/t_{0})^{2}-1|}$ to complete the proof (see
Appendix C for details).

While the proof is constructed around Gaussian random matrices, our model is
flexible enough to accommodate encoding any matrix through the couplings $J$.
This allows for more versatility in proof techniques, that could be lifted in
future work to further strengthen hardness proofs.

###  III.2 Classical hardness of sampling from the distribution

Now, let us discuss the corresponding sampling problem. Following usual
arguments for sampling hardness in Refs. [2, 4, 3], we combine Stockmeyer’s
theorem [44] and an anticoncentration conjecture (see also Appendix F). In
these arguments, Stockmeyer’s theorem provides an
$\mathsf{FBPP}^{\mathsf{NP}}$ algorithm for estimating the output probability
within an inverse polynomial multiplicative error, which is further reduced to
an additive error utilizing the anticoncentration property. Anticoncentration
is characterized by the fact that, for any $\mathbf{x}$, the output
probability $p_{C}(\mathbf{x})$ for a circuit $C$ with $n$ qubits fluctuates
only near $2^{-n}$ over random circuit instances $C$.

However, these techniques cannot be directly applied to our case, as our
problem does not possess the hiding property, unlike those examples. The
hiding property holds when, for a given circuit $C(\theta)$ with a parameter
$\theta$ sampled from a random distribution, the output probability for a
given outcome string $\mathbf{x}$ is the same as that for $\mathbf{0}$ but
with a different parameter, i.e.,
$|\braket{\mathbf{x}}{C(\theta)}{\mathbf{0}}|^{2}=\braket{\mathbf{0}}{C(\theta^{\prime})}{\mathbf{0}}|^{2}$.
As a result, properties of the output distribution averaged over the circuit
parameter $\theta$ do not depend on the outcome bitstring $\mathbf{x}$. For
example, the average-case hardness of
$|\braket{\mathbf{0}}{C(\theta)}{\mathbf{0}}|^{2}$ over $\theta$ holds for
$|\braket{\mathbf{x}}{C(\theta)}{\mathbf{0}}|^{2}$ with any $\mathbf{x}$, and
the anticoncentration property of
$|\braket{\mathbf{0}}{C(\theta)}{\mathbf{0}}|^{2}$ implies that
$|\braket{\mathbf{x}}{C(\theta)}{\mathbf{0}}|^{2}$ also anticoncentrates for
any $\mathbf{x}$.

Compared to those cases, our complexity results (Theorem 1) are only satisfied
for $\mathbf{x}\in X_{m}$ with sufficiently large $m$. The main result
obtained in this section is that such restricted complexity results are
sufficient to prove the difficulty of sampling when the error between the
distribution we sample from and the true distribution $p(\mathbf{x};J;t)$ is
inverse polynomially small, provided that the output distribution
anticoncentrates at least for a constant fraction of $\mathbf{x}\in X_{n/2}$.

Specifically, we obtain the following theorem by applying Stockmeyer’s theorem
[44] to $p(\mathbf{x})$ with $\mathbf{x}\in X_{n/2}$ (see Appendix F for
details).

######  Theorem 2.

Suppose that there exists a classically efficient oracle, $\mathcal{O}$, that
can sample from a probability distribution $p^{\prime}(\mathbf{x})$ which
approximates $p(\mathbf{x})$ within an additive error of $\nu$ in $L_{1}$
norm. We further assume that $p(\mathbf{x})$ and $p^{\prime}(\mathbf{x})$ have
the same domain, i.e., $p^{\prime}(\mathbf{x})=0$ for all $\mathbf{x}\notin X$
when $H=H_{3,4}$. Then for $\gamma\in(0,1)$, there is an
$\mathsf{FBBP}^{\mathsf{NP}}$ machine that estimates $p(\mathbf{x})$ within an
additive error of

  * •

$\epsilon_{I}=(1+g)\gamma^{-1}\nu(2n)2^{-2n}+gp(\mathbf{x})$ ,

  * •

$\epsilon_{II}=(1+g)(\sqrt{\pi}\gamma/2)^{-1}\nu n^{1/2}{2n\choose
n}^{-1}+gp(\mathbf{x})$ ,

for $H_{1,2}$ and $H_{3,4}$, respectively, with a success probability of $\geq
1-\gamma$ for $\mathbf{x}\in X_{n/2}$, where $g=1/\mathrm{poly}(n)$.

We also need the following anticoncentration conjecture for the hardness of
sampling.

######  Conjecture 2 (Anticoncentration).

Let
$p(\mathbf{x};J;t)=|\braket{\mathbf{x}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$
for $\mathbf{y}_{0}$ given in Conjecture 1. Then, there exists
$t_{1}=\Theta(\log n)$ such that for all $t>t_{1}$, the following holds for a
constant fraction, $\zeta$, of $\mathbf{x}\in X_{n/2}$:

  * •

There exist universal constants $0<\alpha,\beta<1$ such that

$\displaystyle\underset{J\sim\mathcal{N}(0,1)^{n\times
n}}{\mathrm{Pr}}\Bigl{[}p(\mathbf{x};J;t)>\alpha 2^{-2n}\Bigr{]}\geq\beta,$
(6)

  * •

There exist universal constants $0<\alpha,\beta<1$ such that

$\displaystyle\underset{J\sim\mathcal{N}(0,1)^{n\times
n}}{\mathrm{Pr}}\Bigl{[}p(\mathbf{x};J;t)>\alpha{2n\choose
n}^{-1}\Bigr{]}\geq\beta,$ (7)

for $H_{1,2}$ and $H_{3,4}$, respectively.

Here, $t_{1}=\Theta(\log n)$ is obtained from the numerical observation
presented in the next section. As we consider $U=e^{-iHt}$ with $\|H\|=O(n)$,
our condition states that anticoncentration holds when $\|H\|t=O(n\log n)$.
This is consistent with the condition that the sparse IQP demonstrates
anticoncentration [45], the circuit of which can be written as
$e^{-iH^{\prime}}$ with $\|H^{\prime}\|=O(n\log n)$.

Mark et al. [46] provides further evidence for Conjecture 2 at least for
$H_{2,3,4}$ and large $t$. Specifically, the study showed that the output
distribution converges to the Porter–Thomas distribution as
$t\rightarrow\infty$ when the eigenvalues of the Hamiltonian follow a
reasonable condition, which is likely to be satisfied for non-integrable
Hamiltonians. As $H_{2,3,4}$ are non-integrable, their results suggest that
the output distribution of $H_{2,3,4}$ anticoncentrates for large $t$ since
the Porter–Thomas distribution implies anticoncentration. Although this result
only applies to sufficiently large $t$, given that anticoncentration is weaker
than the Porter–Thomas distribution, we expect that the output distribution
can anticoncentrate for smaller $t$. In Appendix G, we also prove a result
supporting anticoncentration for $H_{1}$. We particularly show that the time-
averaged value of $p(\mathbf{x};J;t)$ must be $2\times 2^{-2n}$ for all
$\mathbf{x}\in X$ regardless of $J$. Thus, if $t$ is sufficiently large to
equilibrate all $p(\mathbf{x};J;t)$, the output distribution anticoncentrates.

By combining Theorem 2 and Conjecture 2, we obtain the following corollary
(see Appendix F for details).

######  Corollary 1.

Suppose that there exists $t_{0}=\Theta(\log n)$ such that Conjectures 1 and 2
hold. If there exists a classically efficient algorithm that can sample from a
probability distribution which approximates $p_{1,2}(\mathbf{x};J;t_{0})$
$(p_{3,4}(\mathbf{x};J;t_{0}))$ within an additive error of $\varepsilon\zeta
n^{-1}/8$ $(\varepsilon\zeta\sqrt{\pi}n^{-1/2}/8)$, then there exists an
$\mathsf{FBPP}^{\mathsf{NP}}$ algorithm that can approximate
$p(\mathbf{x};J;t_{0})$ within an additive error of $\varepsilon 2^{-2n}$
$(\varepsilon{2n\choose n}^{-1})$ with probability greater than $1-\beta$ over
$J$. As it is a _#P-hard_ problem according to Conjecture 1, this implies that
the PH collapses to the third level.

Corollary 1 implies that sampling from the output distribution within an
inverse polynomially small error is a classically hard problem unless the PH
collapses. We additionally note that the required $L_{1}$ distance can be even
relaxed to a constant for $H_{3,4}$ if a constant fraction of $\mathbf{x}\in
X$ (instead of $\mathbf{x}\in X_{n/2}$) satisfies the anticoncentration
condition (see Appendix F for details).

Figure 3: Left: Output probability $p(\mathbf{x};J;t)$ averaged over
$J\sim\mathcal{N}(0,1)^{n\times n}$ and $\mathbf{x}\in X_{n/2}$ as a function
of $t$. Results are plotted for $n=4$ (the lightest) to $n=10$ (the darkest)
with a step size of $2$. The plot indicates that time to equilibrate $t_{\rm
eq}$ scales with $\log n$. Right: Scatter plots of
$\mathbb{E}_{J}[p(\mathbf{x};J;t)]{2n\choose n}$ (x-axis) versus
$\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]{2n\choose n}^{2}$ (y-axis). For $t/\log
n\in[1,2,3,4]$ and $n\in[4,6,8,10]$, we numerically compute
$p(\mathbf{x};J;t)$ and $p(\mathbf{x};J;t)^{2}$ for all $\mathbf{x}\in
X_{n/2}$ averaged over $J\sim\mathcal{N}(0,1)^{n\times n}$. Each data point
represents $\mathbb{E}_{J}[p(\mathbf{x};J;t)]{2n\choose n}$ and
$\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]{2n\choose n}^{2})$ for each
$\mathbf{x}$. The dashed lines indicate
$\mathbb{E}_{J}[p(\mathbf{x})]{2n\choose n}=1/2$ and
$\mathbb{E}_{J}[p(\mathbf{x})^{2}]{2n\choose n}^{2}=4$. We also present $r$
which is the ratio of $\mathbf{x}\in X_{n/2}$ satisfying
$\mathbb{E}_{J}[p(\mathbf{x})]{2n\choose n}\geq 1/2$ and
$\mathbb{E}_{J}[p(\mathbf{x})^{2}]{2n\choose n}^{2}\leq 4$ (see main text).
Numerical results are obtained by fully diagonalizing the Hamiltonian for
$n\in[4,6,8]$ and by using the second-order Trotter–Suzuki formula with a time
step of $dt=10^{-6}$ for $N=10$. We selected such a small time step to ensure
that the Trotter error is sufficiently smaller than $p(\mathbf{x};J;t)$
itself, which scales as $2^{-2n}$. The results are averaged over $2^{10}$ and
$2^{8}$ random instances of $J\sim\mathcal{N}(0,1)^{n\times n}$ for
$n\in[4,6,8]$ and $n=10$, respectively.

##  IV Numerical verification of anticoncentration

In this section, we provide numerical evidence of the anticoncentration
property given in Conjecture 2. We mainly focus on the XX model ($H_{3}$) here
and provide results for the Ising model ($H_{1}$) in Appendix G.

To numerically test anticoncentration, it is sufficient to compute
$\mathbb{E}_{J}[p(\mathbf{x};J;t)]$ and
$\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]$ for all $\mathbf{x}\in X_{n/2}$. Let
us suppose that $\mathbf{x}$ satisfies the following conditions:

$\displaystyle\mathbb{E}_{J}[p(\mathbf{x};J;t)]$ $\displaystyle\geq
K{2n\choose n}^{-1},$ (8) $\displaystyle\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]$
$\displaystyle\leq\Lambda{2n\choose n}^{-2},$ (9)

where $K\leq 1$ and $\Lambda\geq 1$ are arbitrary constants. Applying Paley-
Zygmund inequality to $p(\mathbf{x};J;t)$ yields the following

$\displaystyle\underset{J\sim\mathcal{N}(0,1)^{n\times
n}}{\mathrm{Pr}}\Bigl{[}p(\mathbf{x};J;t)>\theta\mathbb{E}_{J}[p(\mathbf{x};J;t)]\Bigr{]}$
$\displaystyle\qquad\geq(1-\theta)^{2}\frac{\mathbb{E}_{J}[p(\mathbf{x};J;t)]^{2}}{\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]},$
(10)

which is satisfied for all $0\leq\theta\leq 1$. By entering $\theta=1/2$, we
obtain

$\displaystyle\underset{J\sim\mathcal{N}(0,1)^{n\times
n}}{\mathrm{Pr}}\Bigl{[}p(\mathbf{x};J;t)\geq\frac{K}{2}{2n\choose
n}^{-1}\Bigr{]}\geq\frac{K^{2}}{4\Lambda}.$ (11)

Comparing the above equation to Conjecture 2, we obtain $\alpha=K/2$ and
$\beta=K^{2}/(4\Lambda)$. Therefore, if the ratio of $\mathbf{x}\in X_{n/2}$
satisfying Eq. (9) converges to a constant as $n$ increases, we numerically
confirm Conjecture 2. For numerical results, we choose $K=1/2$ and
$\Lambda=4$, yielding $\alpha=1/4$ and $\beta=1/64$.

For $n\in[4,6,8,10]$, we plot the output probability $p(\mathbf{x};J;t)$
averaged over $\mathbf{x}\in X_{n/2}$ and $J\sim\mathcal{N}(0,1)^{n\times n}$
in Fig. 3 (Left). The results show that the averaged output probability
reaches a stationary value for $t\gtrsim 3\log n$ regardless of $n$. This
resembles a typical behavior of quantum equilibration [47]. However, compared
to most existing literature on the equilibration of local observables, we
observe the equilibration of the output probability $p(\mathbf{x};J;t)$.

We also present $p(\mathbf{x};J;t)$ and $p(\mathbf{x};J;t)^{2}$ averaged over
$J$ in Fig. 3 (Right). For each $n$ and $t$ presented in the plot, we also
compute $r$, the ratio of $\mathbf{x}\in X_{n/2}$ which satisfies Eq. (9). Our
results suggest that $r$ does not converge to zero as $n\rightarrow\infty$ for
any of $t/\log n\in[2,3,4]$. Specifically, we observe that $r$ oscillates with
$n$ for $t/\log n\in[2,3]$, and $r\geq 0.7$ for $t=4\log n$ regardless of $n$.
For any value of considered $t$, $r$ is unlikely to converge to zero with $n$,
thus suggesting the anticoncentration property by the Paley-Zygmund
inequality.

We remark that the anticoncentration results we presented here rule out the
possibility that $p(\mathbf{x};J;t)$ is close to $0$ for most of $\mathbf{x}$.
In this sense, anticoncentration is not merely an ingredient for the hardness
proof but is also regarded as a signature of the average-case hardness to some
extent (see Ref. [1] for further discussion).

##  V Quantum circuit realizations

So far, we demonstrated that classical algorithms for computing the output
distribution or sampling from it do not exist unless the PH collapses. In this
section, we show that a quantum computer can perform the sampling task in
polynomial time using a quantum circuit implementation of the spin dynamics
(i.e., $U=e^{-\mathrm{i}Ht}$).

As some experimental platforms such as trapped ions [48, 49] and neutral atoms
[50, 51] naturally implement spin-$1/2$ Hamiltonians with noisy qubits, it is
tempting to implement our Hamiltonians directly on these systems. However, an
important obstacle is that our Hamiltonians are defined on a bipartite
lattice, which makes a direct implementation challenging without adjusting the
connectivity of the Hamiltonian. While such connectivity-adjusted Hamiltonians
may maintain classical hardness thanks to some degrees of long-range
connectivity in those systems, it is an open question that needs separate
proof. Noisy qubits are another limitation of using such native systems. With
the recent advances in simulating quantum systems using tensor networks,
quantum advantage experiments with noisy qubits, once believed to be
classically intractable, became amenable to classical simulations [52, 53].

Hence, we discuss an implementation based on fault-tolerant intermediate-scale
quantum (ISQ) hardware. Most importantly, we prove that a circuit with a
polynomial depth is sufficient to generate samples with a small $L_{1}$
distance, which is classically challenging by Corollary 1.

As the unitary operator generated by $H_{1}$,
$e^{-\mathrm{i}H_{1}t}=\prod_{ij}e^{-\mathrm{i}(tJ_{ij}/n)\sigma^{(i)}_{x}\tau^{(j)}_{x}}$,
can be directly implemented in a quantum circuit with $n^{2}$ logical gates
due to its commuting nature, we here focus on $H_{2,3,4}$ with non-commuting
terms. Still, we note that there is an interesting relation between the
dynamics generated by $H_{1}$ and IQP circuits [4, 45, 22, 23]. They both can
be implemented using commuting gates, but the underlying hardness problems are
different: the permanent for $H_{1}$ and the universality of measurement-based
quantum computation for IQP. Thus, while we expect that most classical
algorithms for IQP circuits can also simulate the dynamics of $H_{1}$, the
required number of qubits to demonstrate quantum advantage counted based on
the underlying problems can be much smaller for the dynamics of $H_{1}$ (see,
e.g., Ref. [54] and our discussion below).

We now prove that a polynomial number of gates is sufficient to demonstrate
quantum advantage for Hamiltonians $H_{2,3,4}$. As a first step, we show that
implementing a unitary $\tilde{U}$ that approximates $U=e^{-\mathrm{i}Ht_{0}}$
within an inverse polynomially small additive error suffices to demonstrate
quantum advantage. For $\rho=\ket{\mathbf{y}_{0}}\bra{\mathbf{y}_{0}}$, the
output probability distributions from two unitary operators are given as
$p(\mathbf{x})=\braket{\mathbf{x}}{U\rho U^{\dagger}}{\mathbf{x}}$ and
$p^{\prime}(\mathbf{x})=\braket{\mathbf{x}}{\tilde{U}\rho\tilde{U}^{\dagger}}{\mathbf{x}}$.
Then, the $L_{1}$ distance between the two probability distributions is given
by

$\displaystyle\sum_{\mathbf{x}}|p(\mathbf{x})-p^{\prime}(\mathbf{x})|$
$\displaystyle=\sum_{\mathbf{x}}|\braket{\mathbf{x}}{U\rho
U^{\dagger}}{\mathbf{x}}-\braket{\mathbf{x}}{\tilde{U}\rho\tilde{U}^{\dagger}}{\mathbf{x}}|$
$\displaystyle\leq 2\operatorname{Tr}|U\rho
U^{\dagger}-\tilde{U}\rho\tilde{U}^{\dagger}|$ $\displaystyle\leq
2\operatorname{Tr}|U\rho(U^{\dagger}-\tilde{U}^{\dagger})|+2\operatorname{Tr}|(U-\tilde{U})\rho\tilde{U}^{\dagger}|$
$\displaystyle\leq 4\|U-\tilde{U}\|.$

For the first inequality, we used the property of the trace distance
$2\operatorname{Tr}|\rho-\rho^{\prime}|=\max_{\\{\Pi_{i}\\}}\sum_{i}\operatorname{Tr}[\Pi_{i}(\rho-\rho^{\prime})]$,
in which the maximization is over all POVM $\\{\Pi_{i}\\}$. The second and
last inequalities are obtained from the triangular inequality and from
Hölder’s inequality $\operatorname{Tr}|AB|\leq\|A\|\operatorname{Tr}|B|$ with
$\operatorname{Tr}|\rho U|=1$ which is satisfied for any density matrix $\rho$
and any unitary operator $U$, respectively.

We now consider a quantum circuit implementation of
$U=e^{-\mathrm{i}Ht_{0}}=(e^{-\mathrm{i}Ht_{0}/M})^{M}$ by using the Trotter
decomposition [55]. By applying the $p$-th order Suzuki–Trotter formula, we
obtain a quantum circuit implementation of $e^{-\mathrm{i}Ht_{0}/M}$ with an
additive error of $O[\Upsilon_{p}(t/M)^{p+1}]$, where

$\displaystyle\Upsilon_{p}=\sum_{i_{1},\cdots,i_{p+1}=1}^{\Gamma}\|[h_{i_{p+1}},\cdots,[h_{i_{2}},h_{i_{1}}]\|,$
(12)

with $H=\sum_{i=1}^{\Gamma}h_{i}$. As we apply $M$ repetition of
$e^{-\mathrm{i}Ht_{0}/M}$, the total error becomes
$O[M\Upsilon_{p}(t_{0}/M)^{p+1}]$

For example, let us compute the error from the first order ($p=1$) Trotter
decomposition using $H=H_{3}$. From the Suzuki–Trotter formula, we obtain

$\displaystyle\exp\biggl{[}-i\frac{t_{0}}{M}\sum_{ij}\frac{J_{ij}}{2n}\Bigl{(}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{y}\tau^{(j)}_{y}\Bigr{)}\biggr{]}=$
$\displaystyle\prod_{ij}\exp\biggl{[}-i\frac{t_{0}}{M}\frac{J_{ij}}{2n}\Bigl{(}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{y}\tau^{(j)}_{y}\Bigr{)}\biggr{]}+O(\Upsilon_{1}(t_{0}/M)^{-2}),$

where

$\displaystyle\Upsilon_{1}=\sum_{i,j,k,l=1}^{n}\frac{|J_{ij}J_{kl}|}{4n^{2}}\times$
$\displaystyle\qquad\Bigl{\|}\bigl{[}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{y}\tau^{(j)}_{y},\sigma^{(k)}_{x}\tau^{(l)}_{x}+\sigma^{(k)}_{y}\tau^{(l)}_{y}\bigr{]}\Bigr{\|}.$
(13)

The commutator is nonzero only when either $i=k$ or $j=l$ (but not both);
therefore, we have $2(n-1)$ possible choices of $(k,l)$ which contribute to
the sum for each $(i,j)$. As there exists a total of $n^{2}$ possibilities in
selecting $(i,j)$, we have $\Upsilon_{1}=O(n)$. Thus by taking
$M=\Theta(n^{2}\log^{2}n)$, we can bound the error by
$O(\Upsilon_{1}M^{-1}t_{0}^{2})=O(n^{-1})$ for $t_{0}=\Theta(\log n)$, which
is sufficient to show the classical hardness by Corollary 1. In this case, the
total number of gates is $Mn^{2}=\Theta(n^{4}\log^{2}n)$.

Similarly, the same type of calculation for the second-order Trotter formula
results in $\Upsilon_{2}=O(n)$, the total additive error of
$O(nM^{-2}t_{0}^{3})$, and $2Mn^{2}$ total number of gates. By choosing
$M=\Theta(n\log^{3/2}n)$, we can bound the error by $O(n^{-1})$ using
$\Theta(n^{3}\log^{3/2}n)$ gates in total. We can generalize this further by
using the $p$-th order formula with arbitrary $p$; however, in any case, we
can bound the error by $O(n^{-1})$ with a polynomial number of gates.

Finally, we roughly estimate how many logical qubits are required to
demonstrate quantum advantage. Recall that our problem’s hardness relies on
the permanent of submatrices, the sizes of which are determined by the output
bitstrings. As the output bitstrings are typically in $X_{m}$ with $m\approx
n/2$, and computing the permanent of a dense matrix, the size of which is
bigger than $\approx 50\times 50$, is considered classically challenging (see,
e.g., Refs. [56, 57]), we believe that an experiment with around $200$ qubits
(i.e., $n\approx 100$) may be sufficient to demonstrate quantum advantage. In
addition, the number of gates to implement the Hamiltonian evolution for
$n=100$ and $t_{0}=5\times\log(100)\approx 23.03$ is approximately given by
$1.2\times 10^{8}$, $3.8\times 10^{8}$, and $1.2\times 10^{9}$, when the
target Trotter error is $10^{-1}$, $10^{-2}$, and $10^{-3}$, respectively (see
Appendix H for details).

However, unlike BosonSampling, where the output probability distribution is
proportional to the permanent of a matrix, our main theorem only states that
computing the output probability distribution implies computing the permanent;
the converse is not necessarily true. Thus, the output probability
distribution of smaller systems can still be classically intractable for
different reasons. Indeed, 50-100 qubits are already regarded as extremely
challenging for usual classical simulation algorithms for quantum spin
dynamics (e.g., tensor network methods, Krylov subspace method, and time-
dependent variational Monte-Carlo) unless the Hamiltonian has more structures
that can be utilized (e.g., geometrically local connectivity).

##  VI Conclusion and Outlook

In this paper, we have proven the classical hardness of quantum spin dynamics.
We showed that the near-exactly computing the output distribution from the
unitary operator generated by a wide class of many-body spin Hamiltonians is
#P-hard. Those Hamiltonians include the Ising, XX, and Heisenberg models,
defined on a bipartite lattice. With the anticoncentration conjecture,
according to which the output probability for each bitstring fluctuates only
near its averaged value over the random Hamiltonian instances, we also proved
that sampling from the output distribution is classically hard. In contrast,
the same task can be completed using a quantum circuit with a polynomial
number of gates.

We expect that our techniques can be applied to other spin models.
Straightforward examples include the Ising, XX, and Heisenberg models on a
bipartite lattice ($H_{1,2,3,4}$ considered in this study) with local fields
in the $z$-direction (see Appendix B). Adding such local fields does not
change $\braket{\mathbf{x}}{H^{l}}{\mathbf{y}_{0}}$ for all $k\leq n$, and the
same complexity arguments hold.

Our results mainly focus on the hardness of computing the output distribution
and the sampling task. Sampling is an essential process for computing
expectation values of many observables, e.g., the spin-$z$ correlation
functions, using a quantum computer. However, the classical hardness of
estimating these observables must be studied separately; our results do not
exclude the existence of an efficient algorithm for approximating
$\braket{\psi_{0}}{e^{iHt}Me^{-iHt}}{\psi_{0}}$, where $M$ is a local
observable. For example, some quantum circuits, including the IQP [58] and 2D
short-depth circuits [59], allow efficient classical algorithms for estimating
local observables, even though sampling from the output distribution is
computationally hard in general. We do not expect an efficient classical
algorithm for estimating general observables for our Hamiltonians to exist, as
our models are non-commuting and have more complex connectivity. However, we
leave it as an open question, given that a formal proof of hardness for
computing observables will require different techniques.

Another open question is the verification scheme for the sampling task. For a
small system in which the time evolution can be classically computed, we can
use the usual cross-entropy benchmarking [3]. However, in contrast to
instantaneous quantum polynomial circuits or random circuits, whose output
distribution becomes the Porter–Thomas distribution, the distribution from our
Hamiltonian dynamics does not necessarily follow it. Hence, a detailed
verification scheme for a large number of qubits must be analyzed.

We provided a preliminary analysis of the requirements for an experiment
demonstrating quantum advantage. These estimates are encouraging in suggesting
that such experiments may be possible using a few hundred error-corrected
qubits with less than a billion gates; these are within expectations of the
capabilities of the first generation of fault-tolerant quantum computers.
However, we emphasize that a more careful and detailed analysis is required to
establish these numbers more confidently.

We close our paper by commenting on the generalization of our theory for going
beyond the bipartite graph, which is essential for practical applications in
chemistry, biology, and materials science. We expect to draw similar
conclusions for non-bipartite graphs, which generate hafnians or loop-hafnians
instead of permanents. Computing hafnians and loop-hafnians is #P-hard as well
(see, e.g., Ref. [60]), and we can repeat the proving processes as done in
this paper. The hardness results for general graphs would close the gap
between our theory and quantum advantage in practical applications, for
example, nuclear magnetic resonance spectroscopy [61], where the complexity of
estimating observables must be additionally considered.

###### Acknowledgements.

CYP thanks Prof. Soonwon Choi for bringing Ref. [46] to our attention. JH
thanks Minhyeok Kang for insightful discussions. This work was partly
supported by the Basic Science Research Program through the National Research
Foundation of Korea (NRF), funded by the Ministry of Education, Science and
Technology (NRF-2021M3E4A1038308, NRF-2021M3H3A1038085,
NRF-2022M3H3A106307411, NRF-2023M3K5A1094805, and NRF-2023M3K5A109481311) and
Institute for Information & communications Technology Promotion (IITP) grant
funded by the Korea government(MSIP) (No. 2019-0-00003, Research and
Development of Core technologies for Programming, Running, Implementing and
Validating of Fault-Tolerant Quantum Computing System). This research used
resources of the National Energy Research Scientific Computing Center, a DOE
Office of Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award
NERSC DDR-ERCAP0025705. JH acknowledges Xanadu for hosting his sabbatical year
visit. The source code used for the simulations is available in the GitHub
repository: https://github.com/XanaduAI/hardness-spin-dynamics.

## References

  * Hangleiter and Eisert [2023] D. Hangleiter and J. Eisert, Computational advantage of quantum random sampling, Rev. Mod. Phys. 95, 035001 (2023).
  * Aaronson and Arkhipov [2013] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, Theory Comput. 9, 143 (2013).
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Characterizing quantum supremacy in near-term devices, Nat. Phys. 14, 595 (2018).
  * Bremner  _et al._ [2016] M. J. Bremner, A. Montanaro, and D. J. Shepherd, Average-case complexity versus approximate simulation of commuting quantum computations, Phys. Rev. Lett. 117, 080501 (2016).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Quantum computational advantage using photons, Science 370, 1460 (2020).
  * Madsen  _et al._ [2022] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. Bulmer, F. M. Miatto, L. Neuhaus, L. G. Helt, M. J. Collins, _et al._ , Quantum computational advantage with a programmable photonic processor, Nature 606, 75 (2022).
  * Bromley  _et al._ [2020] T. R. Bromley, J. M. Arrazola, S. Jahangiri, J. Izaac, N. Quesada, A. D. Gran, M. Schuld, J. Swinarton, Z. Zabaneh, and N. Killoran, Applications of near-term photonic quantum computers: software and algorithms, Quantum Sci. Technol. 5, 034010 (2020).
  * Bravyi  _et al._ [2022] S. Bravyi, O. Dial, J. M. Gambetta, D. Gil, and Z. Nazario, The future of quantum computing with superconducting qubits, J. Appl. Phys. 132 (2022).
  * Beverland  _et al._ [2022] M. E. Beverland, P. Murali, M. Troyer, K. M. Svore, T. Hoefler, V. Kliuchnikov, G. H. Low, M. Soeken, A. Sundaram, and A. Vaschillo, Assessing requirements to scale to practical quantum advantage, arXiv preprint arXiv:2211.07629 (2022).
  * Fetter and Walecka [1971] A. L. Fetter and J. D. Walecka,  _Quantum theory of many-particle systems_ (McGraw-Hill, 1971).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Physical review 140, A1133 (1965).
  * Ashcroft and Mermin [1976] N. Ashcroft and N. Mermin,  _Solid State Physics_, HRW international editions (Holt, Rinehart and Winston, 1976).
  * Levine [2014] I. Levine,  _Quantum Chemistry_ , Pearson advanced chemistry series (Pearson, 2014).
  * Avella  _et al._ [2012] A. Avella, F. Mancini, _et al._ ,  _Strongly Correlated Systems: Numerical Methods_ (Springer, 2012).
  * Elfving  _et al._ [2020] V. E. Elfving, B. W. Broer, M. Webber, J. Gavartin, M. D. Halls, K. P. Lorton, and A. Bochevarov, How will quantum computers provide an industrially relevant computational advantage in quantum chemistry?, arXiv preprint arXiv:2009.12472 (2020).
  * Kitaev  _et al._ [2002] A. Y. Kitaev, A. Shen, and M. N. Vyalyi,  _Classical and quantum computation_ , Graduate Studies in Mathematics No. 47 (American Mathematical Soc., 2002).
  * Bookatz [2014] A. D. Bookatz, Qma-complete problems, Quantum Info. Comput. 14, 361–383 (2014).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * Gao  _et al._ [2017] X. Gao, S.-T. Wang, and L.-M. Duan, Quantum supremacy for simulating a translation-invariant ising spin model, Phys. Rev. Lett. 118, 040502 (2017).
  * Bermejo-Vega  _et al._ [2018] J. Bermejo-Vega, D. Hangleiter, M. Schwarz, R. Raussendorf, and J. Eisert, Architectures for quantum simulation showing a quantum speedup, Phys. Rev. X 8, 021010 (2018).
  * Haferkamp  _et al._ [2020] J. Haferkamp, D. Hangleiter, A. Bouland, B. Fefferman, J. Eisert, and J. Bermejo-Vega, Closing gaps of a quantum advantage with short-time hamiltonian dynamics, Phys. Rev. Lett. 125, 250501 (2020).
  * Fefferman  _et al._ [2017] B. Fefferman, M. Foss-Feig, and A. V. Gorshkov, Exact sampling hardness of Ising spin models, Phys. Rev. A 96, 032324 (2017).
  * Peropadre  _et al._ [2017] B. Peropadre, A. Aspuru-Guzik, and J. J. García-Ripoll, Equivalence between spin Hamiltonians and boson sampling, Phys. Rev. A 95, 032327 (2017).
  * Clinton  _et al._ [2022] L. Clinton, T. Cubitt, B. Flynn, F. M. Gambetta, J. Klassen, A. Montanaro, S. Piddock, R. A. Santos, and E. Sheridan, Towards near-term quantum simulation of materials, arXiv preprint arXiv:2205.15256 (2022).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Evidence for the utility of quantum computing before fault tolerance, Nature 618, 500 (2023).
  * Tindall  _et al._ [2023] J. Tindall, M. Fishman, M. Stoudenmire, and D. Sels, Efficient tensor network simulation of IBM’s kicked ising experiment, arXiv preprint arXiv:2306.14887 (2023).
  * Kechedzhi  _et al._ [2023] K. Kechedzhi, S. Isakov, S. Mandrà, B. Villalonga, X. Mi, S. Boixo, and V. Smelyanskiy, Effective quantum volume, fidelity and computational cost of noisy quantum processing experiments, arXiv preprint arXiv:2306.15970 (2023).
  * Begušić and Chan [2023] T. Begušić and G. K. Chan, Fast classical simulation of evidence for the utility of quantum computing before fault tolerance, arXiv preprint arXiv:2306.16372 (2023).
  * Bouland  _et al._ [2019] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani, On the complexity and verification of quantum random circuit sampling, Nat. Phys. 15, 159 (2019).
  * Movassagh [2023] R. Movassagh, The hardness of random quantum circuits, Nat. Phys. , 1 (2023).
  * Bouland  _et al._ [2022] A. Bouland, B. Fefferman, Z. Landau, and Y. Liu, Noise and the frontier of quantum supremacy, in  _2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)_ (IEEE, 2022) pp. 1308–1317.
  * Krovi [2022] H. Krovi, Average-case hardness of estimating probabilities of random quantum circuits with a linear scaling in the error exponent, arXiv preprint arXiv:2206.05642 (2022).
  * [35] J. M. Arrazola, From NISQ to ISQ, https://pennylane.ai/blog/2023/06/from-nisq-to-isq/, posted: 06-21-2023.
  * Bourassa  _et al._ [2021] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, _et al._ , Blueprint for a scalable photonic fault-tolerant quantum computer, Quantum 5, 392 (2021).
  * Campbell [2021] E. T. Campbell, Early fault-tolerant simulations of the hubbard model, Quantum Sci. Technol. 7, 015007 (2021).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Practical quantum advantage in quantum simulation, Nature 607, 667 (2022).
  * Liesen and Strakos [2013] J. Liesen and Z. Strakos,  _Krylov subspace methods: principles and analysis_ (Numerical Mathematics and Scie, 2013).
  * Arora and Barak [2009] S. Arora and B. Barak,  _Computational complexity: a modern approach_ (Cambridge University Press, 2009).
  * Sherstov [2013] A. A. Sherstov, Making polynomials robust to noise, Theory of Computing 9, 593 (2013).
  * Kane  _et al._ [2017] D. Kane, S. Karmalkar, and E. Price, Robust polynomial regression up to the information theoretic limit, in  _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ (IEEE, 2017) pp. 391–402.
  * Kondo  _et al._ [2022] Y. Kondo, R. Mori, and R. Movassagh, Quantum supremacy and hardness of estimating output probabilities of quantum circuits, in  _2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)_ (IEEE, 2022) pp. 1296–1307.
  * Stockmeyer [1985] L. J. Stockmeyer, On approximation algorithms for #P, SIAM J. Comput. 14, 849 (1985).
  * Bremner  _et al._ [2017] M. J. Bremner, A. Montanaro, and D. J. Shepherd, Achieving quantum supremacy with sparse and noisy commuting quantum computations, Quantum 1, 8 (2017).
  * Mark  _et al._ [2023] D. K. Mark, J. Choi, A. L. Shaw, M. Endres, and S. Choi, Benchmarking quantum simulators using ergodic quantum dynamics, Phys. Rev. Lett. 131, 110601 (2023).
  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, Equilibration, thermalisation, and the emergence of statistical mechanics in closed quantum systems, Rep. Prog. Phys. 79, 056001 (2016).
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nat. Phys. 8, 277 (2012).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, _et al._ , Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, A rydberg quantum simulator, Nat. Phys. 6, 382 (2010).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, _et al._ , Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms, Nature 595, 233 (2021).
  * Pan  _et al._ [2022] F. Pan, K. Chen, and P. Zhang, Solving the sampling problem of the sycamore quantum circuits, Phys. Rev. Lett. 129, 090502 (2022).
  * Oh  _et al._ [2023] C. Oh, M. Liu, Y. Alexeev, B. Fefferman, and L. Jiang, Tensor network algorithm for simulating experimental gaussian boson sampling, arXiv preprint arXiv:2306.03709 (2023).
  * Dalzell  _et al._ [2020] A. M. Dalzell, A. W. Harrow, D. E. Koh, and R. L. La Placa, How many qubits are needed for quantum computational supremacy?, Quantum 4, 264 (2020).
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of trotter error with commutator scaling, Phys. Rev. X 11, 011020 (2021).
  * Wu  _et al._ [2018] J. Wu, Y. Liu, B. Zhang, X. Jin, Y. Wang, H. Wang, and X. Yang, A benchmark test of boson sampling on Tianhe-2 supercomputer, Nat. Sci. Rev. 5, 715 (2018).
  * Lundow and Markström [2022] P. Lundow and K. Markström, Efficient computation of permanents, with applications to Boson sampling and random matrices, J. Comput. Phys. 455, 110990 (2022).
  * Ni and Van den Nest [2013] X. Ni and M. Van den Nest, Commuting quantum circuits: efficiently classical simulations versus hardness results, Quantum Information & Computation 13, 54 (2013).
  * Bravyi  _et al._ [2021] S. Bravyi, D. Gosset, and R. Movassagh, Classical algorithms for quantum mean values, Nat. Phys. 17, 337 (2021).
  * Björklund  _et al._ [2019] A. Björklund, B. Gupt, and N. Quesada, A faster hafnian formula for complex matrices and its benchmarking on a supercomputer, J. Exp. Algorithmics 24, 11 (2019).
  * Sels  _et al._ [2020] D. Sels, H. Dashti, S. Mora, O. Demler, and E. Demler, Quantum approximate Bayesian computation for NMR model inference, Nat. Mach. Intell. 2, 396 (2020).
  * Huh [2022] J. Huh, A fast quantum algorithm for computing matrix permanent, arXiv preprint arXiv:2205.01328 (2022).
  * Choi [2022] M.-S. Choi,  _A Quantum Computation Workbook_ (Springer Cham, 2022).
  * Minsky [1967] M. L. Minsky,  _Computation: Finite and Infinite Machines_ (Prentice-Hall, Inc., USA, 1967).
  * Welch and Berlekamp [1986] L. Welch and E. Berlekamp, Error correction for algebraic block codes, U.S. Patent 4,633,470 (1986).
  * Paturi [1992] R. Paturi, On the degree of polynomials that approximate symmetric boolean functions (preliminary version), in  _Proceedings of the twenty-fourth annual ACM symposium on Theory of computing_ (1992) pp. 468–474.
  * Rakhmanov [2007] E. A. Rakhmanov, Bounds for polynomials with a unit discrete norm, Ann. Math. 165, 55 (2007).
  * Devroye  _et al._ [2018] L. Devroye, A. Mehrabian, and T. Reddad, The total variation distance between high-dimensional gaussians with the same mean, arXiv preprint arXiv:1810.08693 (2018).
  * Valiant [1979] L. G. Valiant, The complexity of computing the permanent, Theor. Comput. Sci. 8, 189 (1979).
  * Brandao  _et al._ [2016] F. G. Brandao, A. W. Harrow, and M. Horodecki, Local random quantum circuits are approximate polynomial-designs, Commun. Math. Phys. 346, 397 (2016).

##  Appendix A Big O and related notations

In this appendix, we summarize big $O$ and related notations. Let $f(n)$ and
$g(n)$ be functions of $n\in\mathbb{N}$. Then, each notation is defined as
follows:

  * •

$f(n)=o(g(n))$ iff for all $k>0$, there exists $n_{0}$ such that
$|f(n)|<kg(n)$ for all $n>n_{0}$.

  * •

$f(n)=O(g(n))$ iff there exist $k>0,n_{0}$ such that $|f(n)|\leq kg(n)$ for
all $n>n_{0}$.

  * •

$f(n)=\omega(g(n))$ iff for all $k>0$, there exists $n_{0}$ such that
$f(n)>kg(n)$ for all $n>n_{0}$.

  * •

$f(n)=\Omega(g(n))$ iff there exist $k>0,n_{0}$ such that $f(n)\geq kg(n)$ for
all $n>n_{0}$.

  * •

$f(n)=\Theta(g(n))$ iff $f(n)=O(g(n))$ and $f(n)=\Omega(g(n))$.

For example, if $f(n)=1/\log n$, for any given $k>0$, $f(n)<k$ for all
$n>e^{k}$. Therefore, $f(n)=o(1)$.

##  Appendix B Relating permanent of a matrix and moments of Hamiltonian

In this appendix, we prove how the permanent of $J$ and its submatrices are
obtained from moments of Hamiltonians. We will use the results in this
appendix extensively in the next appendix to prove the hardness of the output
distribution. First, we show the relation between the permanent and $H_{3}$.
Recall that $H_{3}$ is defined as

$\displaystyle H_{3}$
$\displaystyle=\sum_{i,j=1}^{n}\frac{J_{ij}}{2n}\big{[}\sigma^{(i)}_{x}\tau^{(j)}_{x}+\sigma^{(i)}_{y}\tau^{(j)}_{y}\big{]}=\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\big{[}\sigma^{(i)}_{+}\tau^{(j)}_{-}+\sigma^{(i)}_{-}\tau^{(j)}_{+}\big{]},$
(14)

where $\sigma^{(i)}_{\pm}=(\sigma^{(i)}_{x}\pm i\sigma^{(i)}_{y})/2$ (and the
same for $\tau^{(j)}_{\pm}$).

Following the notation in the main text, we consider binary vectors
$\mathbf{x}$ and $\mathbf{y}$ of length $2n$, and their substrings
$\mathbf{x}^{\sigma}=x_{1}\cdots x_{n}$, $\mathbf{x}^{\tau}=x_{n+1}\cdots
x_{2n}$ (the same for $\mathbf{y}_{\sigma,\tau}$). Then, set $X_{k}$ is
defined as
$X_{k}=\\{\mathbf{x}\in\mathbb{Z}_{2}^{2n}:\operatorname{wt}(\mathbf{x^{\sigma}})=k\text{
and }\operatorname{wt}(\mathbf{x^{\tau}})=n-k\\}$. Note that there is only a
single vector in $X_{0}$, which is
$\mathbf{y}_{0}=\mathbf{0}^{\sigma}\mathbf{1}^{\tau}$.

We apply $H_{3}$ to $\ket{\mathbf{y}}$ in $\mathbf{y}\in X_{k}$. In
particular, note that each term of $H_{3}$ generates either a product state
(with some abuse of notation, we use a product state and a bitstring
interchangeably) in $X_{k+1}$ or $X_{k-1}$ with
$J_{ij}\sigma^{(i)}_{+}\tau^{(j)}_{-}$ or with
$J_{ij}\sigma^{(i)}_{-}\tau^{(j)}_{+}$, respectively, from $\ket{\mathbf{y}}$.
In other words, we can write

$\displaystyle H_{3}\ket{\mathbf{y}}=\sum_{\mathbf{x}\in
X_{k+1}|\mathbf{y}}a_{\mathbf{x}}\ket{\mathbf{x}}+\sum_{\mathbf{x}\in
X_{k-1}|\mathbf{y}}b_{\mathbf{x}}\ket{\mathbf{x}},$ (15)

where $X_{k\pm 1}|\mathbf{y}$ denotes a bitstring in $X_{k\pm 1}$ generated
form $\mathbf{y}$ (by applying $J_{ij}\sigma^{(i)}_{\pm}\tau^{(j)}_{\mp}$,
respectively). Given that $\mathbf{y}_{0}\in X_{0}$, Eq. (15) implies that

$\displaystyle\braket{\mathbf{x}}{H_{3}^{l}}{\mathbf{y}_{0}}=0,$ (16)

for all $\mathbf{x}\in X_{k}$ when $l<k$.

Now, let us focus on $\braket{\mathbf{x}}{H_{3}^{k}}{\mathbf{y}_{0}}$ for
$\mathbf{x}\in X_{k}$. For convenience, let us define
$\mathcal{S}=\\{i\in[1,\cdots,n]:x_{i}=1\\}$, and
$\mathcal{T}=\\{i\in[1,\cdots,n]:x_{n+i}=0\\}$. Thus, for $\mathbf{x}\in
X_{k}$, $|\mathcal{S}|=|\mathcal{T}|=k$. By explicitly writing
$\braket{\mathbf{x}}{H_{3}^{k}}{\mathbf{y}_{0}}$, we obtain

$\displaystyle\braket{\mathbf{x}}{H_{3}^{k}}{\mathbf{y}_{0}}=\bra{\mathbf{x}}\overbrace{\Bigl{[}\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\big{(}\sigma^{(i)}_{+}\tau^{(j)}_{-}+\sigma^{(i)}_{-}\tau^{(j)}_{+}\big{)}\Bigr{]}\cdots\Bigl{[}\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\big{(}\sigma^{(i)}_{+}\tau^{(j)}_{-}+\sigma^{(i)}_{-}\tau^{(j)}_{+}\big{)}\Bigr{]}}^{k\text{
times }}\ket{\mathbf{y}_{0}}.$ (17)

Among all possible products, we want to extract terms with $\prod_{i\in
S}\prod_{j\in T}\sigma^{(i)}_{+}\tau^{(j)}_{-}$. As we apply $H_{3}$ for
$k=|\mathcal{S}|=|\mathcal{T}|$ times, terms with only
$\sigma^{(i)}_{+}\tau^{(j)}_{-}$ for $i\in\mathcal{S}$ and $j\in\mathcal{T}$
can contribute to this quantity.

In summary, we have

$\displaystyle\braket{\mathbf{x}}{H^{k}}{\mathbf{y}_{0}}=\frac{1}{n^{k}}\sum_{p\in
S_{k}}\sum_{p^{\prime}\in S_{k}}\prod_{i\in p(\mathcal{S})}\prod_{j\in
p^{\prime}(\mathcal{T})}J_{ij}=\frac{k!}{n^{k}}\sum_{p\in
S_{k}}\prod_{i=1}^{k}J_{\mathcal{S}[i],\mathcal{T}[p(i)]}=\frac{k!}{n^{k}}\mathrm{Per}(J_{\mathcal{ST}}).$
(18)

where $S_{k}$ is the symmetric group with $k$ elements and $p(\mathcal{S})$
for $p\in S_{k}$ is a permutation of $\mathcal{S}$ considering $\mathcal{S}$
as an ordered set. In addition, $J_{\mathcal{ST}}$ is a submatrix of $J$ only
remaining rows in $\mathcal{S}$ and columns in $\mathcal{T}$.

We can extend this logic to $H=H_{4}$, which is defined as

$\displaystyle H_{4}$
$\displaystyle=\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\big{[}\sigma^{(i)}_{+}\tau^{(j)}_{-}+\sigma^{(i)}_{-}\tau^{(j)}_{+}+\sigma^{(i)}_{z}\sigma^{(j)}_{z}\big{]}.$
(19)

The main difference to the previous Hamiltonian is that applying $H_{4}$ to
$\ket{\mathbf{y}}$ for $\mathbf{y}\in X_{k}$ also generates terms in $X_{k}$,
i.e.,

$\displaystyle H_{4}\ket{\mathbf{y}}=\sum_{\mathbf{x}\in
X_{k+1}|\mathbf{y}}a_{\mathbf{x}}\ket{\mathbf{x}}+\sum_{\mathbf{x}\in
X_{k-1}|\mathbf{y}}b_{\mathbf{x}}\ket{\mathbf{x}}+\sum_{\mathbf{x}\in
X_{k}|\mathbf{y}}c_{\mathbf{x}}\ket{\mathbf{x}}.$ (20)

However, because the last term does not affect
$\braket{\mathbf{x}}{H_{4}^{l}}{\mathbf{y}_{0}}$ for all $l<k$ when
$\mathbf{x}\in X_{k}$, we still must apply $H_{4}$ at least $k$ times to
obtain non-vanishing terms of $\mathbf{x}\in X_{k}$ from $\mathbf{y}_{0}\in
X_{0}$.

Similarly, we can compute $\braket{\mathbf{x}}{H_{4}^{k}}{\mathbf{y}_{0}}$. By
using the aforementioned argument, only coefficients for
$\sigma^{(i)}_{+}\tau^{(j)}_{-}$ with $i\in\mathcal{S},j\in\mathcal{T}$
contribute to this quantity. Therefore, we also find that

$\displaystyle\braket{\mathbf{x}}{H_{4}^{k}}{\mathbf{y}_{0}}=\frac{k!}{n^{k}}\mathrm{Per}(J_{\mathcal{ST}}).$
(21)

Let us now consider $H_{1}$, which can be reformulated as

$\displaystyle
H_{1}=\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\sigma^{(i)}_{x}\tau^{(j)}_{x}=\sum_{i,j=1}^{n}\frac{J_{ij}}{n}\Bigl{[}\sigma^{(i)}_{+}\tau^{(j)}_{+}+\sigma^{(i)}_{+}\tau^{(j)}_{-}+\sigma^{(i)}_{-}\tau^{(j)}_{+}+\sigma^{(i)}_{-}\tau^{(j)}_{-}\Bigr{]}.$
(22)

Unlike previous Hamiltonians, applying $H_{1}$ to $\mathbf{y}\in X_{m}$ also
creates terms not in $X:=\cup_{i}X_{i}$ but with the Hamming weights of $n+2$
and $n-2$. Nevertheless, the overall results remain the same because these
terms only go back to $X_{m}$ when another $H_{1}$ is applied. In other words,
if we apply $H_{1}$ twice, we obtain chains, such as
$X_{k}\rightarrow(n+2)\rightarrow X_{k}$ or
$X_{k}\rightarrow(n+2)\rightarrow(n+4)$ where $(\cdot)$ indicates the Hamming
weight of the resulting state. Similarly, applying $H_{1}$ by $l$ times only
creates terms in $X_{k-l}$ to $X_{k+l}$ and some other terms outside of $X$,
i.e.,

$\displaystyle H_{1}^{l}\ket{\mathbf{y}}=\sum_{\mathbf{x}\in
X_{k-l}\cup\cdots\cup
X_{k+l}|\mathbf{y}}a_{\mathbf{x}}\ket{\mathbf{x}}+\sum_{\mathbf{x}\notin\cup_{i}X_{i}}b_{\mathbf{x}}\ket{\mathbf{x}},$
(23)

for $\mathbf{y}\in X_{k}$. This implies
$\braket{\mathbf{x}}{H^{l}}{\mathbf{y}_{0}}$ is still $0$ for all $l<k$ if
$\mathbf{x}\in X_{k}$.

Again, for $\braket{\mathbf{x}}{H_{1}^{k}}{\mathbf{y}_{0}}$, only terms with
$\sigma^{(i)}_{+}\tau^{(i)}_{-}$ contribute to the result, which yields

$\displaystyle\braket{\mathbf{x}}{H_{1}^{k}}{\mathbf{y}_{0}}=\frac{k!}{n^{k}}\mathrm{Per}(J_{\mathcal{ST}}),$
(24)

for $\mathbf{x}\in X_{k}$ as earlier examples. (See also Huh [62]; the
relation can be tested easily with symbolic computing packages such as Q3
[63], written in Wolfram language.)

Finally, we also obtain the same result for $H_{2}$, as
$\sigma^{(i)}_{z}\tau^{(j)}_{z}$ terms cannot change any result for
$\braket{\mathbf{x}}{H^{l}}{\mathbf{y}_{0}}$ with $l\leq k$ (when
$\mathbf{x}\in X_{k}$).

To understand how the power of Hamiltonian is related to the permanent more
easily, we provide an example of a puzzle involving two baskets $\sigma$ and
$\tau$, with $n$ red and $n$ blue balls, each labeled from $1$ to $n$. A
constraint is that baskets $\sigma$ and $\tau$ can only contain red and blue
balls, respectively. There are four available operations:

  * •

A: Remove a blue ball $j$ from $\tau$ and add a red ball $i$ to $\sigma$
($\sigma^{(i)}_{+}\tau^{(j)}_{-}$).

  * •

B: Remove a red ball $i$ from $\sigma$ and add a blue ball $j$ to $\tau$
($\sigma^{(i)}_{-}\tau^{(j)}_{+}$).

  * •

C: Remove a red ball $i$ and a blue ball $j$ from each basket $\sigma$ and
$\tau$ ($\sigma^{(i)}_{-}\tau^{(j)}_{-}$).

  * •

D: Add a red ball $i$ and a blue ball $j$ to each basket $\sigma$ and $\tau$
($\sigma^{(i)}_{+}\tau^{(j)}_{+}$).

We start from the initial state that no ball is in $\sigma$ and $n$ balls are
in $\tau$. If we can use exactly $k$ operations to ensure $k$ balls in
$\sigma$ and $n-k$ balls in $\tau$, the only possible way is to apply A by $k$
times. If we further require all balls with label $S$ be in $\sigma$ and all
balls with $T$ be removed from $\tau$, the balls for operation A must be
selected as $i\in\mathcal{S}$ for $\sigma$ and $j\in\mathcal{T}$ for $\tau$.
For each path from the initial state to the final state, the product of
$J_{ij}$ for $(i,j)$ from each application of A contributes to
$\braket{\mathbf{x}}{H_{1}^{k}}{\mathbf{y}_{0}}$. As all permutations of
$\mathcal{S}$ and $\mathcal{T}$ are allowed for each path, their summation is
given by Eq. (18).

We remark that the same argument is also valid even when local fields in the
$z$-direction are present, i.e., for

$\displaystyle\tilde{H}_{1,2,3,4}=H_{1,2,3,4}+\sum_{i}h^{(1)}_{i}\sigma^{i}_{z}+\sum_{j}h^{(2)}_{j}\tau^{(j)}_{z},$
(25)

one can check that

$\displaystyle\braket{\mathbf{x}}{\tilde{H}_{1,2,3,4}^{l}}{\mathbf{y}_{0}}=0\text{
for all
}l<k,\quad\braket{\mathbf{x}}{\tilde{H}_{1,2,3,4}^{k}}{\mathbf{y}_{0}}=\frac{k!}{n^{k}}\mathrm{Per}(J_{\mathcal{ST}}),$
(26)

when $\mathbf{x}\in X_{k}$. Thus, the results in the main text can also be
applied to these Hamiltonians, implying that estimating output probabilities
for these Hamiltonians is also classically hard (when $h_{i}^{(1)}=O(1)$ and
$h_{j}^{(2)}=O(1)$).

##  Appendix C Complexity of the output distribution

This appendix details the proof of Theorem 1 discussed in the main text with
additional exact hardness results. We mainly consider
$\mathbf{x}_{0}=\mathbf{1}^{\sigma}\mathbf{0}^{\tau}\in X_{n}$ for hardness
results, with some comments for general $\mathbf{x}\in X_{m}$ in Appendix C.8.

###  C.1 Polynomial expansion of $p(\mathbf{x}_{0};J;t)$

We prove several hardness results of estimating $p(\mathbf{x}_{0};J;t)$
throughout this section. Our main tool is a polynomial expansion of
$p(t)=|\braket{\mathbf{x}_{0}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$ (where
we exclude $\mathbf{x}_{0}$ and $J$ from the arguments when they are obvious),
which contains the permanent of $J$, i.e.,
$p(t)=\mathrm{Per}(J)^{2}t^{2n}/n^{2n}+\cdots$. Let us define
$f(t)=\braket{\mathbf{x}_{0}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}$ the Taylor
expansion of which is given by

$\displaystyle f(t)=\sum_{k=0}^{\infty}a_{k}(-\mathrm{i}t)^{k},\quad
a_{k}=\frac{\braket{\mathbf{x}_{0}}{H^{k}}{\mathbf{y}_{0}}}{k!}.$ (27)

From the results in Appendix B, we know that
$\braket{\mathbf{x}_{0}}{H^{k}}{\mathbf{y}_{0}}=0$ for all $k<n$ and
$\braket{\mathbf{x}_{0}}{H^{n}}{\mathbf{y}_{0}}=n!\mathrm{Per}(J)/n^{n}$ for
any $H=H_{1,2,3,4}$ Therefore, the polynomial expansion of $p(t)=|f(t)|^{2}$
is given by

$\displaystyle p(t)=\sum_{k=0}^{\infty}a_{k}t^{k},$ (28)

where $a_{k}=0$ for all $k<2n$ and $a_{2n}=\mathrm{Per}(J)^{2}/n^{2n}$. Thus,
computing the coefficient for $t^{2n}$ of $p(t)$ (or equivalently,
$p^{(2n)}(0)/(2n)!$ where $p^{(2n)}(t)$ is the $2n$-th derivative of $p(t)$)
is classically difficult under Conjecture 3 introduced in the next subsection.

###  C.2 Gaussian matrix permanent conjecture

We introduce the Gaussian matrix permanent conjecture, which is central to our
results.

######  Conjecture 3 (Real-valued version of Theorem 7 in Ref. [2]).

For constants $\varepsilon,\delta\in(0,1)$, estimating $\mathrm{Per}(J)^{2}$
within an additive error of $\varepsilon n!$ with a probability of at least
$1-\delta$ over $J\sim\mathcal{N}(0,1)^{n\times n}$ is _#P-Hard_.

Aaronson and Arkhipov [2] provided two main arguments why Conjecture 3 is
plausible. First, they showed that estimating $\mathrm{Per}(J)^{2}$ within an
additive error of $\varepsilon n!$ is polynomial-time equivalent to estimating
$\mathrm{Per}(J)^{2}$ within a multiplicative error of $\varepsilon$ under the
permanent anticoncentration conjecture (PAC). Since estimating the permanent
within a small multiplicative error is widely believed to be #P-hard, and
there are multiple pieces of evidence that the PAC is true, estimating the
permanent within an additive error of $\varepsilon n!$ is also likely to be
#P-hard. Second, they proved that estimating the permanent of Gaussian
matrices within an additive error of $2^{-\mathrm{poly}(n)}$ with probability
$3/4+1/\mathrm{poly}(n)$ is indeed a #P-hard problem. See also Appendix E for
further discussion.

###  C.3 Hardness of exactly computing the output distribution for
exponentially short time

Conjecture 3 readily implies that the exact evaluation of $p(t)$ at $t\leq
n^{-cn}$ with any $c>1/2$ is classically difficult. Herein, we prove this
fact.

Let us recall the polynomial expansion of
$f(t):=\braket{\mathbf{x}_{0}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}$. By
Taylor’s theorem, there exists $0<\tau_{0}<t$ such that

$\displaystyle
f(t)=\frac{\mathrm{Per}(J)}{n^{n}}(-\mathrm{i}t)^{n}+\frac{f^{(n+1)}(\tau_{0})}{(n+1)!}t^{n+1},$
(29)

where

$\displaystyle
f^{(n+1)}(\tau_{0})=(-\mathrm{i})^{n+1}\braket{\mathbf{x}_{0}}{H^{n+1}e^{-\mathrm{i}H\tau_{0}}}{\mathbf{y}_{0}}.$
(30)

Therefore, $p(t)=|f(t)|^{2}$ becomes

$\displaystyle
p(t)=\frac{t^{2n}}{n^{2n}}\Bigl{|}\mathrm{Per}(J)+\xi(t)\Bigr{|}^{2}$ (31)

where

$\displaystyle\xi(t)=-\mathrm{i}\frac{n^{n}\braket{\mathbf{x}_{0}}{H^{n+1}e^{-\mathrm{i}H\tau_{0}}}{\mathbf{y}_{0}}}{(n+1)!}t,$
(32)

which results in

$\displaystyle|\xi(t)|\leq\frac{\|H\|^{n+1}n^{n}}{(n+1)!}t.$ (33)

Thus, the exact computation of $p(t)$ implies that we can compute

$\displaystyle\mathrm{Per}(J)^{2}+2\mathcal{R}[\xi(t)]\mathrm{Per}(J)+|\xi(t)|^{2},$
(34)

where $\mathcal{R}[\xi(t)]$ is the real part of $\xi(t)$. This implies that we
can estimate $\mathrm{Per}(J)^{2}$ within an additive error of $\epsilon n!$
when $|\xi(t)|=o(\sqrt{n!})$ and $|\mathrm{Per}(J)|=O(\sqrt{n!})$.

From Chebyshev’s inequality, we obtain

$\displaystyle\underset{J\sim\mathcal{N}(0,1)^{n\times
n}}{\mathrm{Pr}}\Bigl{[}|\mathrm{Per}(J)|<k\sqrt{n!}\Bigr{]}\geq
1-\frac{1}{k^{2}},$ (35)

given that the variance of $\mathrm{Per}(J)$ for Gaussian random matrices is
$n!$ [2]. Therefore, by choosing, e.g., $k=10$, we can estimate the permanent
within $o(n!)$ with probability $\geq 0.99$ over random Gaussian matrices $J$
when $|\xi(t)|^{2}=o(n!)$, which is a #P-hard problem by Conjecture 3.

By Stirling’s approximation, for $\|H\|\leq 3n$, we obtain

$\displaystyle|\xi(t)|\leq(3e)^{n+1}n^{n}t,$ (36)

which becomes $o(\sqrt{n!})$ when $t=n^{-cn}$ with $c>1/2$. As $\|H\|\leq 3n$
is satisfied for most of $H$ (see Lemma 9 in Appendix D), it is #P-hard to
exactly compute $p(t)$ for small $t$.

We note that a similar result was obtained for the Ising model in Ref. [24].
However, the exact computation considered in this section and Ref. [24] is
impossible for any realistic computational model as a quantity is subject to
an error introduced by finite-precision arithmetic [64]. Moreover, if we
consider sampling tasks (see Sec. III.2), $t$ must be sufficiently large (at
least $\Omega(\log n)$) to make the output distribution anticoncentrate. Thus,
we must consider how the error introduced in the estimation of $p(t)$ for
larger $t$ affects our estimation of $\mathrm{Per}(J)^{2}$.

###  C.4 Average-case hardness of exactly computing the truncated output
distributions

We now consider the average-case hardness with respect to $t$. Based on the
findings in Ref. [31], we consider a classical oracle that exactly computes a
truncated approximation of
$p(t)=|\braket{\mathbf{x}_{0}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$. In
other words, an oracle returns an exact value of $p_{K}(t)$ given as

$\displaystyle
p_{K}(t)=\sum_{k=0}^{K}\frac{p^{(k)}(0)}{k!}t^{k}=\frac{\mathrm{Per}(J)^{2}}{n^{2n}}t^{2n}+\cdots+\frac{p^{(K)}(0)}{K!}t^{K},$
(37)

with $K\geq 2n$. Here, $p^{(K)}(t)$ is the $K$-th derivative of $p(t)$. Then,
we can prove the hardness of exactly computing $p_{K}(t)$ by using the
Berlekamp–Welch algorithm.

######  Theorem 3 (Berlekamp–Welch algorithm [65]).

Let $q(t)$ be a degree $d$ polynomial. Suppose a dataset
$\\{t_{i},y_{i}\\}_{i=1}^{L}$ with all distinct $\\{t_{i}\\}$ with the promise
that $y_{i}=q(t_{i})$ at least $\min(d+1,(L+d)/2)$ points. Then, $q(t)$ can be
exactly recovered in $\mathrm{poly}(L,d)$ deterministic time.

Now, suppose that there exists an oracle that exactly computes $p_{K}(t)$ over
$t\in[a,b]$ for any $b-a>0$ with probability $>1/2+1/\mathrm{poly}(n)$. By the
Chernoff bound, we can find $L=\mathrm{poly}(n)$ such that at least $(L+d)/2$
points are correct with a success probability of $1-e^{-\Theta(n)}$.
Consequently, we can recover $p_{K}(t)$ by the theorem.

Given that we recovered $p_{K}(t)$ exactly, all coefficients of $p_{K}(t)$ are
also known (in fact, the Berlekamp-Welch algorithm returns all the
coefficients). By extracting the coefficient for $t^{2n}$, we can solve the
square of the permanent of $J$, which is a #P-hard problem by Conjecture 3.
This result can be summarized by the following theorem.

######  Theorem 4.

For any $K\geq 2n$, it is _#P-hard_ to exactly compute $p_{K}(t)$ with
probability $1/2+1/\mathrm{poly}(n)$ over choices of $t\in[a,b]$ with $a<b$.

###  C.5 Approximate hardness from robust polynomial regression

Theorem 4 only proves the average-case hardness for an exact oracle returning
$p_{K}(t)$. Given that our function $p(t)$ is not an exact polynomial
function, our goal is to find a polynomial that approximates $p(t)$ assuming
that data points $\\{t_{i},p(t_{i})\\}$ can be obtained from the oracle. An
approach suggested in Refs. [2, 31] is to combine Paturi’s Lemma [66] and
Rakhmanov’s result [67]. In this approach, the error from the estimation of
$p(t)$ obtains an additional factor, which is exponentially large in the
degree of the polynomial. It can be compensated by requiring exponentially
small estimation errors.

Parturi’s lemma and Rakhmanov’s results are presented as follows.

######  Lemma 1 (Parturi [66]).

Let $q(t)$ be a polynomial of degree $d$, and suppose $|q(t)|\leq\delta$ for
all $|t|\leq\Delta$. Then $q(1)\leq\delta\exp(2d(1+\Delta^{-1}))$.

######  Lemma 2 (Simplified version of Rakhmanov [67]).

Let $q(t)$ be a polynomial of degree $d$, and suppose that
$|q(t_{i})|\leq\delta$ for equidistant points $\\{t_{i}\\}_{i=1}^{L}$ in
$[-1,1]$ with $L>d$. Then $|q(t)|\leq O(\epsilon)$ for all $|t|\leq 1$.

By combining those lemmas, it is possible to prove the average-case
approximate hardness of estimating the permanent of Gaussian random matrices
[2] (albeit the same technique does not work for RCS; see Appendix C.7).

This type of technique, which provides error bounds between a polynomial we
want to know (ground truth) and the polynomial reconstructed from noisy
samples, is called robust polynomial regression. By utilizing recent progress
in this field [32, 41, 42, 33, 34, 43], we can prove a similar result for
estimating high-order derivatives as follows.

######  Lemma 3.

Let $q(t)$ be a degree-$d$ polynomial, and suppose that a dataset
$\\{(t_{i},y_{i})\\}_{i=1}^{L}$ is given, in which $\\{t_{i}\\}$ are $L=d+1$
equidistant points in $[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ and
$|y_{i}-q(t_{i})|\leq\delta$ for all $i$, with $0<\Delta<1$. Then, the
coefficient of $q(t)$ for $t^{k}$ can be estimated in $\mathrm{poly}(d)$ time
deterministically within an error of $O[\delta
t_{0}^{-k}(4/\Delta)^{d}{d\choose k}]$.

We prove Lemma 3 in Appendix D by constructing a degree-$d$ polynomial using
the Lagrange interpolation and by bounding the difference between the
constructed polynomial and $q(t)$.

We are ready to prove a primitive version of Theorem 1.

######  Theorem 5 (FP algorithm).

Let $t_{0}=O(\log n)$ and $0<\Delta<1$ be a constant. Under Conjecture 3,
there exists a constant $c_{0}>0$ such that the following statement holds: If
there exists an efficient oracle that evaluates $p(\mathbf{x};J;t)$ for
$\mathbf{x}\in X_{m}$ with $m\geq\sqrt{n}$ and
$t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ within an additive error of
$n^{-c_{0}n}$ with a probability of $\geq 1-1/n^{2}$ over $t$ and a constant
probability of $\eta>0$ over $J\sim\mathcal{N}(0,1)^{n\times n}$, then
$\mathsf{FP}=\mathsf{\\#P}$.

We prove Theorem 5 for $\mathbf{x}=\mathbf{x}_{0}\in X_{n}$ in the rest of
this subsection (see Appendix C.8 for other $\mathbf{x}\in X_{m}$). We again
use $p_{K}(t)$, a truncated Taylor series of $p(t)$ with degree $K$, to
estimate how the error from the estimation of $p(t)$ (where we exclude
$\mathbf{x}_{0}$ and $J$ from the arguments as when they are obvious) affects
our estimation of $p^{(2n)}(0)/(2n)!$. Let us define $p_{K}(t)$ such that

$\displaystyle p_{K}(t)=\sum_{k=2n}^{K}a_{k}t^{k},$ (38)

where $a_{k}=p^{(k)}(0)/k!$ with a special case of
$a_{2n}=\mathrm{Per}(J)^{2}/n^{2n}$. Using Taylor’s theorem, we have

$\displaystyle|p(t)-p_{K}(t)|\leq\frac{|p^{(K+1)}(\tau)|}{(K+1)!}t^{K+1},$
(39)

for $0\leq\tau\leq t$. With $p(t)=f(t)f(t)^{*}$, where
$f(t)=\braket{\mathbf{x}_{0}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}$, we can
expand $p^{(K+1)}(\tau)$ as follows:

$\displaystyle p^{(K+1)}(\tau)=\sum_{l=0}^{K+1}{K+1\choose
l}f^{(l)}(\tau)[f^{(K+1-l)}(\tau)]^{*}.$ (40)

Using an upper bound of $f^{(l)}(\tau)$ as given in below

$\displaystyle|f^{(l)}(\tau)|=|\braket{\mathbf{x_{0}}}{H^{l}e^{-\mathrm{i}H\tau_{0}}}{\mathbf{y}_{0}}|\leq\|H\|^{l},$
(41)

where $\|H\|$ is the operator norm of $H$, we can bound $p^{(K+1)}(\tau)$ as

$\displaystyle|p^{(K+1)}(\tau)|\leq\sum_{l=0}^{K+1}{K+1\choose
l}\|H\|^{K+1}\leq(2\|H\|)^{K+1}.$ (42)

Finally, we have the truncation error bound of probability, denoted
$\epsilon_{K}$ as follows:

$\displaystyle|p(t)-p_{K}(t)|\leq\frac{(2\|H\|t)^{K+1}}{(K+1)!}:=\epsilon_{K}.$
(43)

Thus, by estimating $p(t)$ within an additive error of $\epsilon$, we can also
estimate $p_{K}(t)$ within an additive error of
$\epsilon^{\prime}:=\epsilon+\epsilon_{K}$.

In summary, the problem reduces to the estimation of $2n$-th coefficient of
$p_{K}(t)$ (i.e., for $t^{2n}$), wherein noisy data for $p(t)$ near $t=t_{0}$
are provided. Then, the theorem is proved as follows: (1) For a given oracle
that approximates $p(\mathbf{x}_{0};J;t)$, we call the oracle $K+1$ times for
a given $J$ to generate an equidistant dataset. The oracle successfully
returns all data points with probability $\geq 1-(K+1)/n^{2}$. (2) We use
Lemma 3 to estimate $a_{2n}$, which will be our estimation of
$\mathrm{Per}(J)^{2}/n^{2n}$. Thus, the final error for $\mathrm{Per}(J)^{2}$
is given by $O[(nt_{0}^{-1})^{2n}\epsilon^{\prime}(4/\Delta)^{K}{K\choose
2n}]$, where $n^{2n}$ comes from the denominator of $a_{2n}$, and by
identifying $\delta=\epsilon^{\prime}$, $k=2n$, $d=K$ in Lemma 3. (3) By
bounding this error to $o(n!)$, we can prove that the existence of such an
oracle solves a #P-hard problem by Conjecture 3.

For the proof, we restrict $H$ to have $\|H\|\leq\kappa n$ with $\kappa=3$.
Here, $\kappa$ can be an arbitrary constant that sufficiently reduces the
outlier probability, i.e., $\mathrm{Pr}[\|H\|>\kappa n]=o(1)$. This adds to
the failure probability of our algorithm. In Appendix D, we show that $H$
constructed using a random Gaussian matrix $J$ produces a small outlier
probability given by $e^{-\Theta(n^{2})}$ when $\kappa=3$ for all
$H=H_{1,2,3,4}$.

For those $H$, we show that $K=\alpha n\log n$ estimates $\mathrm{Per}(J)^{2}$
within an additive error of $o(n!)$ when $\epsilon=1/n^{c_{0}n}$ for a
constant $c_{0}$. First, let us bound $\epsilon_{K}$ [Eq. (43)]. As
$t_{0}=O(\log n)$, we can find $C>0$ such that $t_{0}\leq C\log n$ for a
sufficiently large $n$. Then, we obtain

$\displaystyle\epsilon_{K}\leq\frac{[6nt_{0}(1+\Delta)]^{K+1}}{(K+1)!}\leq\frac{(12nt_{0})^{K+1}}{(K+1)!}$
$\displaystyle\quad\leq\Bigl{(}\frac{12ent_{0}}{K+1}\Bigr{)}^{K+1}\leq\Bigl{(}\frac{12eC}{\alpha}\Bigr{)}^{\alpha
n\log n+1},$ (44)

where we use $0<\Delta<1$ for the second inequality and Stirling’s
approximation for the third inequality. In addition, if we choose
$c_{0}\geq-\alpha\log(12eC/\alpha)$ for $\epsilon=n^{-c_{0}n}$, we obtain

$\displaystyle\epsilon^{\prime}=\epsilon+\epsilon_{K}=O\biggl{[}\Bigl{(}\frac{12eC}{\alpha}\Bigr{)}^{\alpha
n\log n}\biggr{]}.$ (45)

Next, we extract the error factor from Lemma 3, by identifying $q(t)=p_{K}(t)$
with $d=K=n\log n$ and $k=2n$. From the entropic bound of the binomial
coefficient given by

$\displaystyle{K\choose 2n}\leq 2^{K\mathcal{H}(2n/K)},$ (46)

where $\mathcal{H}(p)=-p\log_{2}(p)-(1-p)\log_{2}(1-p)$ is the Shannon
entropy, it follows that ${K\choose 2n}\leq 1.5^{K}$ for sufficiently large
$n$. Thus, the error factor from Lemma 3 becomes

$\displaystyle
O\Bigl{[}\epsilon^{\prime}t_{0}^{-2n}\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{K}{K\choose
2n}\Bigr{]}$
$\displaystyle=O\Bigl{[}\epsilon^{\prime}t_{0}^{-2n}\Bigl{(}\frac{6}{\Delta}\Bigr{)}^{\alpha
n\log n}\Bigr{]}=O\biggl{[}\Bigl{(}\frac{72eC}{\alpha\Delta}\Bigr{)}^{\alpha
n\log n}\biggr{]}.$ (47)

By multiplying the remaining factor $n^{2n}$ from the denominator of $a_{2n}$
with the last expression, we obtain the overall error for
$\mathrm{Per}(J)^{2}$ by $O(n^{n[2+\alpha\log(72eC/\alpha\Delta)]})$. This can
be further bounded by $O(n^{\mu n})=o(n!)$ with $\mu<1$ as we can find
$\alpha$ such that $\alpha\log(72eC/\alpha\Delta)<-1$. Our result is
summarized as follows: Suppose there exists a classically efficient oracle
($\in\mathsf{FP}$) that can estimate $p(t)$ within an additive error of
$n^{-c_{0}n}$ for $t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ with a constant
probability $\eta>0$ over $J\sim\mathcal{N}(0,1)^{n\times n}$. Then, we can
estimate $\mathrm{Per}(J)^{2}$ within an additive error of $o(n!)$ by calling
the oracle $\mathrm{poly}(n)$ times. The success probability of the algorithm
over $J\sim\mathcal{N}(0,1)^{n\times n}$ is given by
$\eta-e^{-\Theta(n^{2})}-(\alpha n\log n+1)/n^{2}\geq 7\eta/8$ for
sufficiently large $n$. Then, it follows that FP=#P by Conjecture 3, which
also implies P=NP [40].

However, Theorem 5 required an oracle that can estimate $p(\mathbf{x};J;t)$
for a range of $t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$. Thus, such an oracle
is much stronger than we want for Conjecture 1, which estimates
$p(\mathbf{x};J;t)$ only at $t=t_{0}$. Fortunately, it is also possible to
prove the hardness using an oracle for $t=t_{0}$ if we require smaller
additive errors in the estimation and use a more powerful algorithm.

###  C.6 Proof of Theorem 1 for $\mathbf{x}\in X_{n}$

In this subsection, we prove Theorem 1 for $\mathbf{x}=\mathbf{x}_{0}\in
X_{n}$ by providing an FBPP algorithm that estimates $\mathrm{Per}(J)^{2}$
within a sufficiently small error when an oracle that evaluates
$p(\mathbf{x}_{0};J;t)$ at $t=t_{0}$ is provided. One of the main differences
of the algorithm presented in this subsection compared to the one introduced
in Theorem 5 is that the oracle only requires to return the estimated values
of $p(t)$ at $t=t_{0}$ instead of that of the range
$[t_{0}(1-\Delta),t_{0}(1+\Delta)]$. This is because we can estimate $p(t)$
near $t_{0}$ if the estimation of $p(t_{0})$ is possible for Gaussian random
matrices $J$, which follows from the lemma below.

######  Lemma 4.

An oracle that estimates $p(J;t_{0})$ within an additive error of $\epsilon$
with probability $\eta$ over $J\sim\mathcal{N}(0,1)^{n\times n}$ can estimate
$p(J;t)$ within the same error with a probability of at least
$\eta-(3/2)n\sqrt{|(t/t_{0})^{2}-1|}$.

###### Proof.

Let $\mathcal{O}(J)$ be an oracle that satisfies

$\displaystyle\mathrm{Pr}_{J\sim\mathcal{N}(0,1)}[|\mathcal{O}(J)-p(J;t_{0})|\leq\epsilon]>\eta.$
(48)

Let $E$ be the events where $|\mathcal{O}(J)-p(J;t_{0})|\leq\epsilon$. As the
distribution of $p(J;t)$ for $J\sim\mathcal{N}(0,1)$ is the same as that of
$p(J;t_{0})$ with $J\sim\mathcal{N}(0,(t/t_{0})^{2})$, the difference between
the probability of having $E$ is upper bounded by the total variational
distance, i.e.,

$\displaystyle|\mathrm{Pr}_{J\sim\mathcal{N}(0,(t/t_{0})^{2})}[E]-\mathrm{Pr}_{J\sim\mathcal{N}(0,1)}[E]|\leq\delta,$
(49)

where $\delta$ is the total variational distance between two Gaussian
distributions. This completes the proof as
$\delta\leq(3/2)n\sqrt{|(t/t_{0})^{2}-1|}$ [68]. ∎

The lemma implies that, when $\eta\geq 3/4$ is a constant, by choosing
$\Delta\leq\beta/n^{2}$ with $\beta=1/64$, the oracle can estimate $p(t)$
within an error $\epsilon$ with probability $>1/2$ over
$t\in[t_{0}(1-\Delta),t_{0}(1-\Delta)]$. Then, we can apply the following
lemma to estimate the $2n$-th coefficient of $\tilde{q}(t)$, i.e., $a_{2n}$.

######  Lemma 5.

For a degree-$d$ polynomial $q(t)$, suppose that there exists an efficient
oracle that estimates $q(t)$ within error $\delta$ with a constant probability
larger than $1/2$ over $t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$. Then there
exists an _FBPP_ algorithm that returns the coefficient of $q(t)$ for $t^{k}$,
within an additive error of $O[\delta t_{0}^{-k}(4/\Delta)^{d}{d\choose k}]$
with a success probability of at least $2/3$.

See Appendix D for the proof. We now prove Theorem 1 for
$\mathbf{x}=\mathbf{x}_{0}$, which we recall as follows.

See 1

We follow the same steps for the proof for Theorem 5 to find an error bound
for $\mathrm{Per}(J)^{2}$. Let us bound $\epsilon_{K}$ for $K=\alpha n^{3}\log
n$ and $t_{0}\leq C\log n$, which is given by

$\displaystyle\epsilon_{K}\leq\Bigl{(}\frac{12ent_{0}}{K+1}\Bigr{)}^{K+1}\leq\Bigl{(}\frac{12ent_{0}}{K}\Bigr{)}^{K+1}$
$\displaystyle=\Bigl{(}\frac{12eC}{\alpha
n^{2}}\Bigr{)}\Bigl{(}\frac{12eC}{\alpha}\Bigr{)}^{\alpha n^{3}\log
n}n^{-2\alpha n^{3}\log n}$
$\displaystyle\leq\Bigl{(}\frac{12eC}{\alpha}\Bigr{)}^{\alpha n^{3}\log
n}n^{-2\alpha n^{3}\log n},$ (50)

where the last inequality is satisfied for $n\geq(12eC/\alpha)^{1/2}$. Thus
for any $c_{1}>2\alpha$ and $\epsilon=n^{-c_{1}n^{3}\log n}$, we have
$\epsilon=o(\epsilon_{K})$, i.e.,

$\displaystyle\epsilon^{\prime}=\epsilon+\epsilon_{K}=O\biggl{[}\Bigl{(}\frac{12et_{0}}{\alpha}\Bigr{)}^{\alpha
n^{3}\log n}n^{-2\alpha n^{3}\log n}\biggr{]}.$ (51)

We now apply Lemma 5 by putting $\delta=\epsilon^{\prime}$,
$\Delta\leq\beta/n^{2}$, $d=K$, and $k=2n$. This yields the following additive
error bound:

$\displaystyle\epsilon^{\prime}t_{0}^{-2n}\left(\frac{4n^{2}}{\beta}\right)^{K}{K\choose
2n}\leq\epsilon^{\prime}t_{0}^{-2n}\left(\frac{6}{\beta}\right)^{\alpha
n^{3}\log n}n^{2\alpha n^{3}\log n}.$ (52)

Consequently, the permanent can be estimated within an additive error of

$\displaystyle
O\Biggl{[}\Bigl{(}\frac{n}{t_{0}}\Bigr{)}^{2n}\Bigl{(}\frac{72et_{0}}{\alpha\beta}\Bigr{)}^{\alpha
n^{3}\log n}\Biggr{]},$ (53)

which is $o(n!)$ for $\alpha>(72et_{0})/\beta$. Hence, the existence of an
efficient classical oracle ($\in\mathsf{FP}$) for estimating $p(t_{0})$ within
an additive error of $n^{-c_{1}n^{3}\log n}$ with probability $\geq 3/4$ over
$J\sim\mathcal{N}(0,1)^{n\times n}$ implies that there exists an FBPP
algorithm for estimating the permanent of Gaussian matrices; this also implies
that FBPP=#P.

###  C.7 Difficulty of combining Lemmas 3 and 4

Given Lemma 4, one might ask whether we can combine Lemmas 3 and 4 to obtain a
deterministic algorithm (i.e., not FBPP but FP) for estimating the $k$-th
coefficient of $q(t)$ only using an oracle for $p(\mathbf{x};J;t_{0})$.
However, the main problem of combining these results is that Lemma 3 requires
$L>d$ data points a priori. They are only accessible when the oracle can
output $p(t)$ with a probability greater than $1-1/L$ over
$t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$, and this is the main reason why we
cannot easily combine Lemmas 3 and 4.

To understand the situation better, let us recall the setting for Theorem 5.
The error from the polynomial approximation (the difference between $p_{K}(t)$
and $p(t)$) was given by $\epsilon_{K}=(12net_{0}/K)^{K}$ which is magnified
by a factor of $t_{0}^{-2n}(4/\Delta)^{K}2^{n\mathcal{H}(2n/K)}\leq
t_{0}^{-2n}(6/\Delta)^{K}$ for the derivatives. Thus, the overall error
becomes

$\displaystyle t_{0}^{-2n}\Bigl{(}\frac{72net_{0}}{K\Delta}\Bigr{)}^{K}.$ (54)

Now let us assume that there exists an oracle that computes $p(t_{0})$ with
probability $\eta>0$ over $J\sim\mathcal{N}(0,1)^{n\times n}$. By Lemma 4, we
can estimate $p(t)$ for $t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ with
probability $\eta-(3/2)n\sqrt{[(1+\Delta)^{2}-1]}$. However, the condition
requiring $L$ data points in advance is only satisfied when the oracle outputs
a correct estimation of $p(t)$ for a given $t_{i}$ with a probability greater
than $1-O(1/K)$. To satisfy this condition, we need $\eta=1-O(1/K)$ and
$\Delta=O(1/(n^{2}K^{2}))$. By entering this $\Delta$ into Eq. (54), we have

$\displaystyle O[t_{0}^{K-2n}(72en^{3}K)^{K}],$ (55)

which diverges with $n$ for any $K\geq 2n$.

This problem was the reason why the first average-case hardness proof of the
RCS [31] did not work. The problem is solved in later studies by using a
fixed-degree rational polynomial [32] or using an algorithm within a higher
level of the PH [43, 33, 34]. In the first approach, the author interpolated
the worst-case and average-case instances of random circuits using a fixed-
degree rational polynomial. However, our problem is unlikely to admit a
similar solution as we are unsure whether such an interpolation involving the
permanent of a matrix even exists. Our approach presented in the main text
using Lemma 5 resembles the second approach.

###  C.8 Proof of Theorem 1 for other $\mathbf{x}\in X_{m}$

In previous subsections, we only have considered the output probability
$p(\mathbf{x};J;t)=|\braket{\mathbf{x}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$
for $\mathbf{x}=\mathbf{x}_{0}=\mathbf{1}^{\sigma}\mathbf{0}^{\tau}\in X_{n}$.
In this case, the $2n$-th derivative of $p(\mathbf{x};J;t)$ with respect to
$t$ at $t=0$ is given by the permanent of the full matrix $J$. This subsection
briefly shows that Theorem 1 still holds for other output bitstrings. From
Appendix B, we know that

$\displaystyle\braket{\mathbf{x}}{H^{l}}{\mathbf{y}_{0}}=0\text{ for all
}l<k,\quad\braket{\mathbf{x}}{H^{m}}{\mathbf{y}_{0}}=\frac{1}{n^{m}}\mathrm{Per}(J_{\mathcal{ST}}),$
(56)

for $\mathbf{x}\in X_{m}$, where $J_{\mathcal{ST}}$ is $m\times m$ submatrix
of $J$, the rows and columns of which are obtained from
$\mathcal{S}=\\{i\in[1,\cdots,n]:x_{i}=1\\}$ and
$\mathcal{T}=\\{i\in[1,\cdots,n]:x_{n+i}=0\\}$, respectively. In addition, we
require $m=n^{\Omega(1)}$ (e.g., $n^{0.001}$, $\sqrt{n}$, $n/2$, etc.) to make
the permanent problem #P-hard w.r.t. $n$.

Again, we define $p_{K}(t)$ as a truncated Taylor series of
$p(t)=|\braket{\mathbf{x}}{e^{-\mathrm{i}Ht}}{\mathbf{y}_{0}}|^{2}$ (where we
consider it for fixed $J$ and $\mathbf{x}$), which is given by

$\displaystyle p_{K}(t)=\sum_{k=m}^{K}a_{k}t^{k},\quad
a_{2m}=\frac{1}{n^{2m}}\mathrm{Per}(J_{\mathcal{ST}})^{2}.$ (57)

Then, the difference between $p_{K}(t)$ and true polynomial $p(t)$ is given by

$\displaystyle|p(t)-p_{K}(t)|\leq\frac{(2\|H\|t)^{K+1}}{(K+1)!}:=\epsilon_{K},$
(58)

which is the same as before, i.e., the oracle estimates $p_{K}(t)$ within an
additive error of $\epsilon+\epsilon_{K}$. Now, we reconstruct
$\hat{p}(t)=\sum_{k=2m}^{K}\hat{a}_{k}t^{k}$ by using a dataset obtained from
the oracle that returns an estimated value of $p(t)$, and let our estimation
of $\mathrm{Per}(J_{\mathcal{ST}})^{2}$ be $n^{2m}\hat{a}_{2m}$. Collecting
all the error factors, our estimation error from the recovery algorithms is
written as:

$\displaystyle\Bigl{(}\frac{n}{t_{0}}\Bigr{)}^{2m}(\epsilon+\epsilon_{K})\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{K}{K\choose
2m}.$ (59)

Further, according to the Gaussian permanent conjecture, we know that the
estimation of $\mathrm{Per}(J_{\mathcal{ST}})^{2}$ within an additive error of
$\varepsilon m!$ is #P-hard when $m=\mathrm{poly}(n)$. Thus, we aim to find a
condition that the error given in Eq. (59) becomes smaller than $\varepsilon
m!$.

Let us first consider the FP algorithm (Theorem 5). Following the proof for
$\mathbf{x}=\mathbf{x}_{0}$, we use $t_{0}\leq C\log n$ for a constant $C$,
$K=\alpha n\log n$, and $t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ with
$\Delta<1$. Then, we obtain

$\displaystyle\epsilon_{K}\leq\Bigl{(}\frac{12eC}{\alpha}\Bigr{)}n^{\alpha
n\log(12eC/\alpha)}.$ (60)

By choosing $c_{0}>\alpha\log(12eC/\alpha)$ from $\epsilon=n^{-c_{0}n}$, we
have $\epsilon_{K}\leq\epsilon$ for sufficiently large $n$, i.e.,
$\epsilon+\epsilon_{K}\leq O[n^{-\alpha n}\log(12eC/\alpha)]$. Next, we also
have

$\displaystyle{K\choose 2m}\leq{K\choose 2n}\leq 1.5^{K},$ (61)

for sufficiently large $n$. By combining all these factors, we have

$\displaystyle\Bigl{(}\frac{n}{t_{0}}\Bigr{)}^{2m}(\epsilon+\epsilon_{K})\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{K}{K\choose
2m}=O\biggl{[}n^{2m}\Bigl{(}\frac{72eC}{\alpha\Delta}\Bigr{)}^{\alpha n\log
n}\biggr{]}.$ (62)

For any $0\leq m\leq n$, we can find $\alpha$ such that
$\alpha\log(72eC/\alpha\Delta)<-1$, which makes the bound Eq. (62) be $o(m!)$.
This implies that $\mathrm{Per}(J_{\mathcal{ST}})^{2}$ is estimated within an
additive error of $o(m!)$, which is #P-hard for $m=\mathrm{poly}(n)$.

Next, we follow similar steps for the FBPP algorithm (Theorem 1). We use the
same parameters for the proof of Theorem 1, i.e., $\Delta\leq\beta/n^{2}$, and
$K=\alpha n^{3}\log n$ for a constant $\alpha,\beta>0$. We then obtain

$\displaystyle(\epsilon+\epsilon_{K})t_{0}^{-2m}\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{K}{K\choose
2m}\leq
O\Bigl{[}t_{0}^{-2m}\Bigl{(}\frac{72et_{0}}{\alpha\beta}\Bigr{)}^{\alpha
n^{3}\log n}\Bigr{]}.$ (63)

Thus, the overall error will be

$\displaystyle
O\Bigl{[}\Bigl{(}\frac{n}{t_{0}}\Bigr{)}^{2m}\Bigl{(}\frac{72et_{0}}{\alpha\beta}\Bigr{)}^{\alpha
n^{3}\log n}\Bigr{]},$ (64)

which is $o(1)$ for any $\alpha>(72et_{0})/\beta$.

###  C.9 Remarks on Theorem 1

We provide the following remarks for Theorem 1.

  1. 1.

Our proof techniques for Theorem 1 can also be directly used to prove a weaker
version of Conjecture 3. In particular, one can prove that under the BPP
reduction (i.e., such an oracle implying FBPP=#P) that computing
$\mathrm{Per}(J)$ or $\mathrm{Per}(J)^{2}$ with an additive error of
$n^{-O(n)}$ with a probability of $3/4$ over random Gaussian matrices $J$, is
#P-hard. Details are provided in Appendix E.

  2. 2.

Our proof for Theorem 1 is robust against the additive error for the permanent
given in Conjecture 3. Both the FP (Theorem 5) and FBPP (Theorem 1) algorithms
can be feasible even for the $n^{-O(n)}$ (instead of $\varepsilon n!$)
additive error estimation of the permanent by adjusting some constants. Thus,
under the BPP reduction, we can directly use the results in Appendix E, i.e.,
we do not require Conjecture 3.

  3. 3.

An open question is whether the required error bounds for oracles in Theorem 1
can be relaxed. We believe that proving (even a restricted version of)
Conjecture 1 by relaxing error bounds from $n^{-\alpha n^{3}\log n}$ for
Theorem 1 to $\varepsilon 2^{-2n}$ could be impactful. However, this is
challenging under our current setup based on a truncated Taylor series unless
significant improvements are made to the robust polynomial regression
techniques. This is because all related approaches for RCS results in a
similar error bound [43, 33, 34], and the technique we used in Lemma 5 (which
was introduced in Ref. [42]) achieves the information theoretic bound.

  4. 4.

An alternative approach for relaxing the required error bound would involve
the search for a closed-form expression of $p(\mathbf{x};J;t)$ written in a
combination of a finite number of elementary functions involving the permanent
of $J$ or some modified matrices. Such an expression would facilitate
extracting $\mathrm{Per}(J)^{2}$ more robustly, thus improving the error bound
of the permanent estimations. We determined the possibility of writing some
high-order terms from the Taylor expansion as a product of the permanent of
$J$ and other complex functions of $J_{ij}$. Nevertheless, as other high-order
terms do not allow simple decompositions, detailed analyses in this direction
can be considered in future work.

##  Appendix D Proofs of Lemmas used for Theorem 1

###  D.1 Proof of Lemma 3

We prove Lemma 3, which provides an FP algorithm for estimating the
derivatives of a polynomial function when noisy data are provided.

######  Theorem 6 (Restatement of Lemma 3).

Let us assume that $q(t)$ is a degree-$d$ polynomial and $0<\Delta<1$. Suppose
that a dataset $\\{(t_{i},y_{i})\\}_{i=1}^{L}$ such that $\\{t_{i}\\}$ are
$L=d+1$ equidistant points in $[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ with
$\Delta<1$ and $|y_{i}-q(t)|\leq\delta$ for all $i$ is given. Then the $k$-th
coefficient of $q(t)$, i.e., $q^{(k)}(0)/k!$, can be estimated in
$\mathrm{poly}(d)$ time within an error $O[\delta
t_{0}^{-k}(4/\Delta)^{d}{d\choose k}]$.

###### Proof.

We construct a degree-$d$ polynomial, $\hat{q}(t)$, by using the Lagrange
interpolation from the given dataset and select its coefficient for $t^{k}$ as
our estimation of the $k$-th coefficient of $q(t)$, i.e., $q^{(k)}(0)/k!$.
Thus, we aim to bound the error between the $k$-th coefficient of two
polynomials $\hat{q}(t)$ and ${q}(t)$. By writing
$r(t):=q(t)-\hat{q}(t)=\sum_{k=0}^{d}b_{k}t^{k}$, the error we want to bound
is $b_{k}$.

Let us define

$\displaystyle R(s)$ $\displaystyle:=r(t=t_{0}(\Delta
s+1)\bigr{)}=\sum_{k=0}^{d}c_{k}s^{k}.$ (65)

Then, condition $|q(t_{i})-y_{i}|\leq\delta$ for $L$ equidistant points
$t_{i}\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$ implies that $|R(s)|\leq\delta$
for $L$ equidistant points between $s\in[-1,1]$. By using Rakhmanov’s theorem
[67], we obtain

$\displaystyle|R(s)|\leq C\delta,$ (66)

for all $|s|<1$ where $C>0$ is a constant. According to the lemma below, we
have

$\displaystyle\sum_{k}|c_{k}|\leq C4^{d}\delta.$ (67)

We now obtain a relation between $b_{k}$ and $c_{k}$ by inserting
$s=(t/t_{0}-1)/\Delta$ to $R(t)$, which yields

$\displaystyle b_{k}=\frac{1}{t_{0}^{k}}\sum_{l\geq
k}^{d}(-1)^{(l-k)}\frac{c_{l}}{\Delta^{l}}{l\choose k}.$ (68)

Therefore,

$\displaystyle|b_{k}|$
$\displaystyle\leq\frac{1}{t_{0}^{k}}\Bigl{(}\sum_{l\geq
k}^{d}|c_{l}|\Bigr{)}\Delta^{-d}{d\choose k}\leq
C\delta\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{d}t_{0}^{-k}{d\choose k},$ (69)

which completes the proof. ∎

######  Remark.

The proof can be simplified without invoking Rakhmanov’s theorem but by
directly using the setting in Ref. [41]. However, this does not change the
overall scaling behavior.

######  Lemma 6 (Lemma 4.1 in Ref. [41]).

Let $Q(x)=\sum_{i=0}^{d}a_{i}x^{i}$ be a polynomial. Then

$\displaystyle\sum_{k=0}^{d}|a_{k}|\leq 4^{d}\max_{x\in[-1,1]}|Q(x)|.$ (70)

One subtlety of the algorithm is that it returns $\hat{q}(t)$ in a
representation involving Lagrange basis polynomials, i.e., we have
$\hat{q}(t)=\sum_{i=1}^{K+1}\alpha_{i}l_{i}(t)$ where

$\displaystyle l_{i}(t)=\prod_{j\neq i}\frac{t-t_{j}}{t_{i}-t_{j}}.$ (71)

Thus, to extract the $k$-th coefficient of $\hat{q}(t)$, the $k$-th
coefficient of $l_{i}(t)$ must be efficiently obtained. The following Lemma
shows that all the coefficients can be obtained in polynomial time.

######  Lemma 7.

For a given degree-$d$ polynomial, which is represented using its roots
$\\{r_{k}\\}$, i.e., $q(t)=\prod_{k=1}^{d}(t-r_{k})$, one can obtain all the
coefficients in $O(d\log^{2}d)$ time.

###### Proof.

Let us define $q_{1}(t)=\prod_{k=1}^{d/2}(t-r_{k})$ and
$q_{2}(t)=\prod_{k=d/2+1}^{d}(t-r_{k})$. When coefficients of $q_{1}(t)$ and
$q_{2}(t)$ are known, coefficients of $q(t)=q_{1}(t)q_{2}(t)$ can be obtained
using the fast Fourier transformation in $O[d\log(d)]$ time. Then, the
recurrence relation is given as follows:

$\displaystyle T(d)=2T(d/2)+O(d\log(d)),$ (72)

where $T(d)$ is the time complexity for a degree-$d$ polynomial. Solving the
equation gives $T(d)=O[d\log^{2}(d)]$. ∎

As mentioned earlier, one of the limitations of Theorem 6 is that the
algorithm requires a dataset with $L=d+1$ equidistant points, which is only
possible when the oracle correctly outputs with a probability greater than
$1-1/L$. We also note that reducing the required number of samples is possible
using a dataset optimally chosen for Chebyshev polynomials as in Ref. [43].
However, such a technique does not improve the overall scaling of the error.

###  D.2 Proof of Lemma 5

######  Theorem 7 (Restatement of Lemma 5).

Let us assume that $q(t)$ is a degree-$d$ polynomial and $0<\Delta<1$. Suppose
that we can obtain samples of $q(t)$ within error $\delta$ with a constant
probability greater than $1/2$ over uniformly random
$t\in[t_{0}(1-\Delta),t_{0}(1+\Delta)]$, i.e.,

$\displaystyle\mathrm{Pr}[|y_{i}-q(t_{i})|\leq\delta]>1/2,\text{ where
}t_{i}\sim\mathcal{U}_{[t_{0}(1-\Delta),t_{0}(1+\Delta)]}.$ (73)

Then, there exists an _FBPP_ algorithm that returns the $k$-th coefficient of
$q(t)$, i.e., $q^{(k)}(0)/k!$, within an additive error of $O[\delta
t_{0}^{-k}(4/\Delta)^{d}{d\choose k}]$ with a success probability of at least
$2/3$.

We prove the theorem by using the following lemma.

######  Lemma 8 (Ref. [42]).

Let $Q(s)$ be a degree-$d$ polynomial defined in $s\in[-1,1]$ and suppose that
we can obtain samples of $Q(s)$ within error $\delta$ with probability
$\rho>1/2$ where each $s_{i}$ is drawn from the Chebyshev distribution
$D_{c}(s)\sim\frac{1}{\sqrt{1-s^{2}}}$. Then, by using
$O[\frac{d}{\epsilon}\log\frac{d}{\epsilon\chi}]$ number of samples, we can
find a polynomial $\hat{Q}(s)$ that satisfies

$\displaystyle\max_{t\in[-1,1]}|Q(s)-\hat{Q}(s)|\leq(2+\epsilon)\delta,$ (74)

with probability $1-\chi$, where $\epsilon\in(0,1/2)$ is a constant.

On the other hand, if $s_{i}$ is drawn from the uniform distribution instead,
then $O(\frac{d^{2}}{\epsilon^{2}}\log\chi^{-1})$ samples suffice for the same
result. In both cases, the runtime of this algorithm is polynomial in the
sample size.

We prove the theorem by applying this lemma to $Q(s):=q(t_{0}(1+s\Delta))$. By
putting $\epsilon=1/4$ into Lemma 8, we can find $\hat{Q}(s)$ such that

$\displaystyle\max_{x\in[-1,1]}|Q(s)-\hat{Q}(s)|\leq\frac{9}{4}\delta,$ (75)

with a probability greater than $2/3$ by using $O(d^{2})$ samples. Then, we
take the $k$-th coefficient of $\hat{q}(t):=\hat{Q}((t/t_{0}-1)/\Delta)$ as
our estimation of $q^{(k)}(0)/k!$.

By writing $R(s):=Q(s)-\hat{Q}(s)=\sum_{k=0}^{d}c_{k}s^{k}$, we obtain

$\displaystyle
r(t):=q(t)-\hat{q}(t)=R\Bigl{(}\frac{t/t_{0}-1}{\Delta}\Bigr{)}=\sum_{k=0}^{d}b_{k}t^{k}.$
(76)

Again, from

$\displaystyle b_{k}=\frac{1}{t_{0}^{k}}\sum_{l\geq
k}^{d}(-1)^{(l-k)}\frac{c_{l}}{\Delta^{l}}{l\choose k},$ (77)

and Lemma 6, we know that

$\displaystyle|b_{k}|\leq\frac{9\delta}{4t_{0}^{k}}\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{d}{d\choose
k},$ (78)

which is the error from our estimation.

###  D.3 Probability to have $\|H\|\leq 3n$

By taking a term-wise norm, the infinite norm of $H_{1,2,3,4}$ is bounded by

$\displaystyle\|H_{1,3}\|\leq\sum_{ij}|J_{ij}|/n,\quad\|H_{2}\|\leq
2\sum_{ij}|J_{ij}|/n,\quad\|H_{4}\|\leq\frac{3}{2}\sum_{ij}|J_{ij}|/n.$ (79)

Hence, the question is when $|J_{ij}|\leq 3n^{2}$ (for $H_{1,3}$) , $3n^{2}/2$
(for $H_{2}$), and $2n^{2}$ (for $H_{4}$) are satisfied. As each $J_{ij}$ is a
random variable, we can use the lemma below to have

$\displaystyle\mathrm{Pr}[\sum_{ij}|J_{ij}|\leq cn^{2}]\geq
1-2^{n^{2}}e^{-c^{2}n^{2}/2}.$ (80)

By choosing $c=3$ for $H_{1,3}$, $c=3/2$ for $H_{2}$, and $c=2$ for $H_{4}$,
we obtain

$\displaystyle\mathrm{Pr}[\|H_{1,3}\|\leq 3n]\geq
1-\Bigl{(}\frac{2}{e^{9}}\Bigr{)}^{n^{2}},$ (81)
$\displaystyle\mathrm{Pr}[\|H_{2}\|\leq 3n]\geq
1-\Bigl{(}\frac{2}{e^{9/4}}\Bigr{)}^{n^{2}},$ (82)
$\displaystyle\mathrm{Pr}[\|H_{4}\|\leq 3n]\geq
1-\Bigl{(}\frac{2}{e^{4}}\Bigr{)}^{n^{2}}.$ (83)

######  Lemma 9.

Let $\\{X_{i}\\}$ be $n$ independent random Gaussian variables, i.e.,
$X_{i}\sim\mathcal{N}(0,1)$ for all $i\in\\{1,\cdots,n\\}$. Then the following
bound is satisfied:

$\displaystyle\mathrm{Pr}\Bigl{[}\frac{1}{n}\sum_{i=1}^{n}|X_{i}|\geq
t\Bigr{]}\leq 2^{n}e^{-nt^{2}/2}.$ (84)

###### Proof.

Let $X\sim\mathcal{N}(0,1)$. Then the standard Gaussian integration gives the
moment generating function $\mathbb{E}[e^{tX}]=e^{t^{2}/2}$ for
$t\in\mathbb{R}$. In addition,

$\displaystyle\mathbb{E}[e^{t|x|}]\leq\mathbb{E}[e^{t|x|}+e^{-t|x|}]=2e^{t^{2}/2}.$
(85)

We now consider $n$ Gaussian random variables . For any $\lambda>0$, we have

$\displaystyle\mathrm{Pr}\Bigl{[}\frac{1}{n}\sum_{i=1}^{n}|X_{i}|\geq
t\Bigr{]}$
$\displaystyle=\mathrm{Pr}\Biggl{[}\exp\Bigl{(}\frac{\lambda}{n}\sum_{i=1}^{n}|X_{i}|\Bigr{)}\geq
e^{\lambda t}\Biggr{]}$ $\displaystyle\leq e^{-\lambda
t}\mathbb{E}\Bigg{[}\exp\Bigl{(}\frac{\lambda}{n}\sum_{i=1}^{n}|X_{i}|\Bigr{)}\Bigg{]}=e^{-\lambda
t}\prod_{i=1}^{n}\mathbb{E}\Bigg{[}\exp\Bigl{(}\frac{\lambda}{n}|X_{i}|\Bigr{)}\Bigg{]}$
$\displaystyle\leq e^{-\lambda
t}2e^{\lambda^{2}/(2n^{2})}=2^{n}e^{\lambda^{2}/(2n)-\lambda t},$ (86)

where we used the standard Chernoff bound for the first inequality.

By taking $\lambda=nt$, which minimizes the argument of the exponential, we
obtain the desired inequality. ∎

##  Appendix E Average-case approximate hardness of a Gaussian matrix
permanent

The arXiv version of Ref. [2] proved that estimating the permanent of Gaussian
matrices within an exponentially small additive error is difficult (Theorem 62
therein), which supports Conjecture 3. This appendix provides an alternative
proof giving an exact required error bound for estimation.

We consider the following well-known fact that computing permanent of 0/1
matrices is #P-hard for worst-case instances.

######  Theorem 8 (Ref. [69]).

There exists a matrix $X\in\mathbb{R}^{n\times n}$, the elements of which are
all $0$ or $1$ such that computing $\mathrm{Per}(X)$ is _#P-hard_.

Note that this readily implies that computing the permanent of 0/1 matrix
within an additive error of $o(1)$ is also #P-hard given that the result must
be an integer. We now prove that estimating the permanent of a Gaussian matrix
is also hard on average.

######  Theorem 9.

If there exists an efficient classical algorithm that estimates
$\mathrm{Per}(J)$ within an additive error of $n^{-2n(1+\delta)}$ for any
$\delta>0$ with probability $3/4$ over $J\sim\mathcal{N}(0,1)^{n\times n}$,
then _FBPP=#P_.

Following Aaronson and Arkhipov [2], we choose a matrix $X(t):=tX+(1-t)Y$
where $X$ is a 0/1 matrix for the worst-case instance, and $Y$ is a matrix
drawn from a random Gaussian distribution. Based on the definition of the
matrix permanent, $q(t)=\mathrm{Per}(X(t))$ is a degree-$n$ polynomial. Under
this setting, the following lemmas will be used to prove the theorem.

######  Lemma 10.

The total variational distance between $X(t):=tX+(1-t)Y$ for a random Gaussian
matrix $Y$ and $X=(x_{ij})$ is bounded above by $n(3\sqrt{t}+t)$.

###### Proof.

For a given $t$, we can consider
$X(t)\sim\prod_{ij}\mathcal{N}(tx_{ij},(1-t)^{2}):=D_{t}$. Based on results in
Ref. [68], for $D_{t}^{\prime}:=\prod_{ij}\mathcal{N}(tx_{ij},1^{2})$, we have

$\displaystyle\|D_{t}-D_{t}^{\prime}\|\leq
n\frac{3}{2}\sqrt{1-(1-t)^{2}}<3n\sqrt{t}.$ (87)

Note that Ref. [2] used a different bound given by
$\|D_{t}-D_{t}^{\prime}\|\leq n^{2}t$ (Lemma 48, therein). However, as we
could not follow that proof, we here use a weaker bound proved in Ref. [68],
which does not change the overall scaling. In addition, we have

$\displaystyle\|D_{t}^{\prime}-\mathcal{N}(0,1)^{n\times
n}\|\leq\sqrt{\sum_{ij}(tx_{ij})^{2}}\leq nt,$ (88)

as shown in Ref. [2]. Then,

$\displaystyle\|D_{t}-\mathcal{N}(0,1)^{n\times n}\|$
$\displaystyle\leq\|D_{t}-D_{t}^{\prime}\|+\|D_{t}^{\prime}-\mathcal{N}(0,1)^{n\times
n}\|$ $\displaystyle\leq n(3\sqrt{t}+t).$ (89)

∎

######  Lemma 11 (Ref. [34]).

Let us assume that $q(t)$ is a degree-$d$ polynomial and $\Delta\in(0,1)$.
Suppose that we can obtain samples of $q(t)$ within an additive error of
$\delta$ with a constant probability of $\eta>1/2$ for $t_{i}$ drawn uniformly
from $[-\Delta,\Delta]$. Then, there exists an _FBPP_ algorithm that returns
$q(1)$ within an additive error of $\frac{9\delta}{4}(4/\Delta)^{d}$.

###### Proof.

Let us define $Q(s):=q(\Delta s)$, a degree-$d$ polynomial in $s$. Under the
condition used in the lemma, an oracle can estimate $Q(s)$ for $s\in[-1,1]$
within an additive error of $\delta$. By using such an oracle, we can
construct a polynomial $\hat{Q}(s)$ following Lemma 8, which satisfies

$\displaystyle\max_{s\in[-1,1]}|Q(s)-\hat{Q}(s)|\leq\frac{9}{4}\delta.$ (90)

Then for $\hat{q}(t)=\hat{Q}(t/\Delta)$, our estimation of $q(1)$ will be
$\hat{q}(1)$.

We now bound the error between $q(1)$ and $\hat{q}(1)$. Let us define
$R(s):=Q(s)-\hat{Q}(s)=\sum_{k=0}^{d}c_{k}s^{k}$. By Lemma 6, we obtain

$\displaystyle\sum_{k=0}^{d}|c_{k}|\leq\frac{9\delta}{4}4^{d}.$ (91)

Then, it follows that

$\displaystyle|q(1)-\hat{q}(1)|=|R(\Delta^{-1})|\leq\sum_{k=0}^{d}|c_{k}|\Delta^{-k}\leq\frac{9\delta}{4}\Bigl{(}\frac{4}{\Delta}\Bigr{)}^{d}.$
(92)

∎

Proof of Theorem 9 Let us assume an oracle that estimates $\mathrm{Per}(J)$
within an additive error of $\epsilon$, with a constant probability of $3/4$
over a random Gaussian matrix $J$. Then from Lemma 10, such an oracle also
estimates $q(t)=\mathrm{Per}(X(t))$ with probability
$\eta=3/4-n(3\sqrt{\Delta}+\Delta)$ when $t\in[-\Delta,\Delta]$. By taking
$\Delta=(16n)^{-2}$, we could obtain $\eta>1/2$ for all $n$.

Next, we apply Lemma 11 to $q(t)=\mathrm{Per}(X(t))$. As the degree of $q(t)$
is $n$, and the error is $\delta=\epsilon$, the error for
$q(1)=\mathrm{Per}(J)$ becomes

$\displaystyle\frac{9\epsilon}{4}\Bigl{(}\frac{4}{(16n)^{-2}}\Bigr{)}^{n}=\frac{9\epsilon}{4}(2^{10}n^{2})^{n}.$
(93)

The error is upper bounded by $o(1)$ when $\epsilon=n^{-2n(1+\delta)}$ for any
$\delta>0$. This implies that estimating $q(t)$ at $t=0$ within an additive
error of $\epsilon$ is #P-hard under the BPP reduction.

The same argument as above also applies to $\mathrm{Per}(J)^{2}$ instead of
$\mathrm{Per}(J)$. From Theorem 28 in Ref. [2], we know that the worst-case
instance $X\in\mathbb{R}^{n\times n}$ exists such that computing
$\mathrm{Per}(X)^{2}$ is #P-hard. Then, we consider a polynomial
$X(t):=tX+(1-t)Y$ where $Y$ is a random matrix from $\mathcal{N}(0,1)^{n\times
n}$. As $\mathrm{Per}(X(t))^{2}$ is a degree-$2n$ polynomial, we can use the
same steps to obtain the following theorem:

######  Theorem 10.

If there exists an efficient classical algorithm that estimates
$\mathrm{Per}(J)^{2}$ within an additive error of $n^{-4n(1+\delta)}$ for any
$\delta>0$ with probability $3/4$ over $J\sim\mathcal{N}(0,1)^{n\times n}$,
then _FBPP=#P_.

Compared to Theorem 62 in Ref. [2], the required probability of our theorem is
$3/4$ (in fact, any constant probability larger than $1/2$ can be used)
instead of $3/4+1/\mathrm{poly}(n)$ and the required additive error
$n^{-4n(1+\delta)}$ does not depend on the probability.

##  Appendix F Proof details of classical hardness of the sampling task

In Appendix C, we proved that computing the output distribution is a
classically hard problem because it belongs to class #P-hard. In this section,
we further prove that sampling from the output distribution is also a
classically challenging problem unless the PH collapses.

We follow the usual arguments in Refs. [2, 4, 3], which combine Stockmeyer’s
theorem and the anticoncentration property of the output distribution to argue
that sampling from the distribution is also classically hard (unless the PH
collapses). In these arguments, Stockmeyer’s theorem provides
$\mathsf{FBPP}^{\mathsf{NP}}$ algorithm for estimating the output probability
within an inverse polynomial multiplicative error, which is further reduced to
an additive error utilizing the anticoncentration property.

However, as we mentioned in the main text, these techniques cannot be directly
applied to our case, as our complexity results (Theorem 1) are only satisfied
for $\mathbf{x}\in X_{m}$ with sufficiently large $m$. This is related to the
fact that our problem lacks the hiding property. The main result obtained in
this section is that such restricted complexity results are sufficient to
prove the difficulty of sampling when the error between the distribution we
sample from and the true distribution $p(\mathbf{x};J;t)$ is inverse
polynomially small, provided that the output distribution anticoncentrates at
least for a constant fraction of $\mathbf{x}\in X_{n/2}$.

Let us start with Stockmeyer’s approximate counting theorem.

######  Theorem 11 (Stockmeyer [44]).

Suppose that there exists an oracle $\mathcal{O}$ that generates a sample from
the probability distribution $p(\mathbf{x})$ by taking a random bitstring
$\mathbf{r}$, i.e.,

$\displaystyle
p(\mathbf{x})=\underset{\mathbf{r}\in\\{0,1\\}^{\mathrm{poly}(n)}}{\mathrm{Pr}}\bigl{[}\mathcal{O}(\mathbf{r})=\mathbf{x}\bigr{]}.$
(94)

There exists an $\mathsf{FBPP}^{\mathsf{NP}^{\mathcal{O}}}$ machine that, for
a given $\mathbf{x}\in\\{0,1\\}^{n}$, approximates $p(\mathbf{x})$ within a
multiplicative factor of $g=1/\mathrm{poly}(n)$. In other words, for a given
$\mathbf{x}$, the machine outputs $\hat{p}(\mathbf{x})$ such that

$\displaystyle(1-g)\hat{p}(\mathbf{x})\leq
p(\mathbf{x})\leq(1+g)\hat{p}(\mathbf{x}).$ (95)

Then, the following theorem provides an argument that a classical approximate
sampler can estimate $p(\mathbf{x})$ for $\mathbf{x}\in X_{n/2}$.

We now prove Theorem 2 introduced in the main text. See 2

###### Proof.

We first prove the case for $H_{1,2}$. From the definition of $L_{1}$ distance

$\displaystyle\nu=\sum_{\mathbf{x}\in\mathbb{Z}_{2}^{2n}}|p(\mathbf{x})-p^{\prime}(\mathbf{x})|,$
(96)

it follows from Markov’s inequality that

$\displaystyle\underset{\mathbf{x}\in\mathbb{Z}_{2}^{2n}}{\mathrm{Pr}}\Bigl{[}|p(\mathbf{x})-p^{\prime}(\mathbf{x})|\geq\frac{\nu}{\delta}2^{-2n}\Bigr{]}\leq\delta.$
(97)

Next, the probability to have $\mathbf{x}\in X_{n/2}$ among all
$\mathbf{x}\in\mathbb{Z}_{2}^{2n}$ is given by

$\displaystyle\frac{{n\choose n/2}^{2}}{2^{2n}}\geq\frac{1}{2n},$ (98)

where we have used the inequality

$\displaystyle{n\choose k}\geq\sqrt{\frac{n}{8k(n-k)}}2^{n\mathcal{H}(k/n)}.$
(99)

Then, by taking $\delta=\gamma(2n)^{-1}$, more than $1-\gamma$ fraction of
$\mathbf{x}\in X_{n/2}$ satisfies

$\displaystyle|p(\mathbf{x})-p^{\prime}(\mathbf{x})|\leq\gamma^{-1}\nu(2n)2^{-2n}.$
(100)

Theorem 11 implies that if there exists an efficient classical sampler for
$p^{\prime}(\mathbf{x})$, then there exists an
$\mathsf{FBPP}^{\mathsf{NP}^{\mathcal{O}}}$ algorithm that outputs
$\hat{p}(\mathbf{x})$ such that

$\displaystyle|\hat{p}(\mathbf{x})-p^{\prime}(\mathbf{x})|\leq
gp^{\prime}(\mathbf{x}),$ (101)

where $g=1/\mathrm{poly}(n)$. Then the distance between $\hat{p}(\mathbf{x})$
and $p(\mathbf{x})$ for $\mathbf{x}$ satisfying Eq. (100) is given by

$\displaystyle|\hat{p}(\mathbf{x})-p(\mathbf{x})|$
$\displaystyle\leq|\hat{p}(\mathbf{x})-p^{\prime}(\mathbf{x})|+|p(\mathbf{x})-p^{\prime}(\mathbf{x})|$
$\displaystyle\leq
gp^{\prime}(\mathbf{x})+|p(\mathbf{x})-p^{\prime}(\mathbf{x})|$
$\displaystyle\leq(1+g)|p(\mathbf{x})-p^{\prime}(\mathbf{x})|+gp(\mathbf{x})$
$\displaystyle\leq(1+g)\gamma^{-1}\nu(2n)2^{-2n}+gp(\mathbf{x}).$ (102)

Proof for $H_{3,4}$ follows the same step as above, besides that we only
consider $\mathbf{x}\in X$ instead of $\mathbb{Z}_{2}^{2n}$. Recall that
$X=\cup_{i=0}^{n}X_{i}$ is the set of all bitstrings that spans the eigenvalue
zero subspace of $J_{z}$. From the definition of $L_{1}$ distance

$\displaystyle\nu=\sum_{\mathbf{x}\in
X}|p(\mathbf{x})-p^{\prime}(\mathbf{x})|,$ (103)

and Markov’s inequality, we obtain

$\displaystyle\underset{\mathbf{x}\in
X}{\mathrm{Pr}}\Bigl{[}|p(\mathbf{x})-p^{\prime}(\mathbf{x})|\geq\frac{\nu}{\delta}{2n\choose
n}^{-1}\Bigr{]}\leq\delta.$ (104)

Here, ${2n\choose n}=|X|$. In addition, the probability to have $\mathbf{x}\in
X_{n/2}$ among all allowed $x$ (within the same $J_{z}$ eigenspace) is given
by

$\displaystyle\frac{{n\choose n/2}^{2}}{{2n\choose
n}}\geq\frac{\sqrt{\pi}}{2}n^{-1/2},$ (105)

where we have used the inequality

$\displaystyle\sqrt{\frac{n}{8k(n-k)}}2^{n\mathcal{H}(k/n)}\leq{n\choose
k}\leq\sqrt{\frac{n}{2\pi k(n-k)}}2^{n\mathcal{H}(k/n)}.$ (106)

Thus by taking $\delta=\gamma(\sqrt{\pi}/2)n^{-1/2}$, more than $1-\gamma$
fraction of $\mathbf{x}\in X_{n/2}$ satisfies

$\displaystyle|p(\mathbf{x})-p^{\prime}(\mathbf{x})|\leq(\sqrt{\pi}\gamma/2)^{-1}\nu
n^{1/2}{2n\choose n}^{-1}.$ (107)

Using the same steps for $H_{1,2}$, the desired error $\epsilon_{II}$ is
obtained. ∎

One can see that the error $\epsilon_{I}$ and $\epsilon_{II}$ are smaller than
that stated in Conjecture 1 denoted by $\varepsilon 2^{-2n}$ and
$\varepsilon{2n\choose n}^{-1}$, respectively, when $p(\mathbf{x})$ also
scales with the inverse of the dimension, i.e., $2^{-2n}$ and ${2n\choose
n}^{-1}$. A version of this condition, so-called the anticoncentration
property [2, 4], is widely known for the sampling problems and can be
rigorously proven for the IQP [4, 22] and the RCS [70]. We also expect that
this property is satisfied for Hamiltonian dynamics, at least for some
$\mathbf{x}\in X_{n/2}$, which we conjecture as follows:

See 2

When the anticoncentration conjecture holds, we can apply Theorem 2 by
choosing $\gamma=\zeta/2$, which makes at least $\zeta/2$ fraction of
$\mathbf{x}\in X_{n/2}$ must have both properties. Then for such $\mathbf{x}$,
by taking $\nu=\varepsilon\gamma/(4n)$ for $H_{1,2}$ (or
$\nu=\varepsilon(\sqrt{\pi}\gamma/4)n^{-1/2}$ for $H_{3,4}$), one can estimate
$p(\mathbf{x};J)$ within an additive error of
$[\varepsilon(1+g)/2+g\alpha]D^{-1}$ with a probability greater than $1-\beta$
over $J$. Here, $D$ is the dimension of the Hilbert space given by $2^{2n}$
and ${2n\choose n}$ for $H_{1,2}$ and $H_{3,4}$, respectively. In addition, as
$g=1/\mathrm{poly}(n)$, $\varepsilon(1+g)/2+g\alpha\leq\varepsilon$ for
sufficiently large $n$.

In summary, we obtained the following result presented in the main text. See 1

Note that Corollary 1 is robust in the sense that one can change the required
fraction of $\mathbf{x}$ showing the anticoncentration property by adjusting
the $L_{1}$ distance $\nu$. For example, the proof for $H_{3,4}$ is feasible
with $\nu=O(\varepsilon)$ when a constant fraction of $\mathbf{x}\in
X=\cup_{i=0}^{n}X_{i}$ satisfies the anticoncentration condition. On the other
hand, if only $\xi=O(1/n)$ fraction of $\mathbf{x}\in X_{n/2}$ shows the
anticoncentration property, the required $L_{1}$ distance becomes
$\nu=O(\varepsilon n^{-2})$ and $\nu=O(\varepsilon n^{-3/2})$ for $H_{1,2}$
and $H_{3,4}$, respectively.

Figure 4: Left: Output probability $p(\mathbf{x};J;t)$ averaged over
$J\sim\mathcal{N}(0,1)^{n\times n}$ and $\mathbf{x}\in X_{n/2}$ as a function
of $t$. We plot results for $n=4$ (the lightest) to $n=10$ (the darkest) with
step size $2$. The plot indicates that time to equilibrate $t_{\rm eq}$ scales
with $\log n$. Right: Scatter plots of
$2^{2n}\mathbb{E}_{J}[p(\mathbf{x};J;t)]$ (x-axis) versus
$2^{4n}\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]$ (y-axis) for each $\mathbf{x}\in
X_{n/2}$. For $t/\log n\in[1,2,3,4]$ and $n\in[4,6,8,10]$, we numerically
compute $p(\mathbf{x};J;t)$ and $p(\mathbf{x};J;t)^{2}$ for all $\mathbf{x}\in
X_{n/2}$ averaged over $2^{10}$ random instances of
$J\sim\mathcal{N}(0,1)^{n\times n}$. The dashed lines indicate
$2^{2n}\mathbb{E}_{J}[p(\mathbf{x})]=1/2$ and
$2^{4n}\mathbb{E}_{J}[p(\mathbf{x})^{2}]=16$. The plot also presents $r$,
which is the ratio of $\mathbf{x}\in X_{n/2}$ satisfying
$2^{2n}\mathbb{E}_{J}[p(\mathbf{x})]\geq 1/2$ and
$2^{4n}\mathbb{E}_{J}[p(\mathbf{x})^{2}]\leq 16$. For numerical simulations,
we used the expansion
$e^{-\mathrm{i}H_{1}t}=\prod_{i,j=1}^{n}[e^{-i(J_{ij}/n)\sigma^{(i)}_{x}\tau^{(j)}_{x}t}]$.

##  Appendix G Anticoncentration of the Ising model

In the main text, we mainly discussed the anticoncentration of the output
distribution for $H_{3}$. In this appendix, we study the anticoncentration
property for the Ising model ($H_{1}$).

###  G.1 Numerical demonstration of anticoncentration

As discussed in the main text, we demonstrated anticoncentration of the output
distribution for $H_{3}$ by computing $\mathbb{E}_{J}[p(\mathbf{x};J;t)]$ and
$\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]$ where
$p(\mathbf{x};J;t)=|\braket{\mathbf{x}}{e^{-\mathrm{i}H_{1}t}}{\mathbf{y}_{0}}|^{2}$.
Similarly, we also compute the same quantity for $H=H_{1}$, as plotted in Fig.
4.

In Fig. 4(Left), we observe that $\mathbb{E}_{\mathbf{x}\in
X_{n/2}}\mathbb{E}_{J}[p(\mathbf{x};J;t)]$ converges to $\approx 2.0\times
2^{-2n}$ for large $t$ regardless of $n$. We analytically prove in the next
subsection that this result is related to the Hamming weight of $\mathbf{x}$,
which is $n$ for $\mathbf{x}\in X$. In addition, the plot shows that the time
for equilibration scales with $\log n$.

Furthermore, we plot $\mathbb{E}_{J}[p(\mathbf{x};J;t)]$ and
$\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]$ for each $\mathbf{x}\in X_{n/2}$ in
Fig. 4(Right). We also present $r$ which is the ratio of $\mathbf{x}$
satisfying

$\displaystyle\mathbb{E}_{J}[p(\mathbf{x};J;t)]\geq\frac{1}{2}2^{-2n}\land\mathbb{E}_{J}[p(\mathbf{x};J;t)^{2}]\leq
16\times 2^{-4n}.$ (108)

The results show that $r$ converges to a non-zero value as $n$ increases for
all $t/\log n\in[1,2,3,4]$. Then it follows from the Paley-Zygmund inequality
that the output probability anticoncentrates.

###  G.2 Converged output probability of the Ising model

In this subsection, we analytically derive the behavior of the output
probability $p(\mathbf{x};J;t)$ when $t\rightarrow\infty$. Our analytic
expression shows that the converged value of the output probability depends on
the Hamming weight of $\mathbf{x}$. For convenience, we divide $\mathbf{x}$
into $\mathbf{x}^{\sigma}=\\{x_{1},\cdots,x_{n}\\}$ and
$\mathbf{x}^{\tau}=\\{x_{n+1},\cdots,x_{2n}\\}$.

As the first step, we expand
$f(\mathbf{x};J;t):=\braket{\mathbf{x}^{\sigma}\mathbf{x}^{\tau}}{e^{-iHt}}{\mathbf{y}_{0}}$,
which is given by

$\displaystyle f(\mathbf{x};J;t)$
$\displaystyle=\braket{\mathbf{x}^{\sigma}\mathbf{x}^{\tau}}{e^{-iHt}}{\mathbf{y}_{0}}=\braket{\mathbf{x}^{\sigma}\mathbf{x}^{\tau}}{e^{-i\sum_{ij}J_{ij}\sigma^{(i)}_{x}\tau^{(j)}_{x}t}}{\mathbf{y}_{0}}$
$\displaystyle=\braket{\mathbf{x}^{\sigma}\mathbf{x}^{\tau}}{H^{\otimes
2n}e^{-i\sum_{ij}J_{ij}\sigma^{(i)}_{z}\tau^{(j)}_{z}t/n}H^{\otimes
2n}}{\mathbf{y}_{0}}$
$\displaystyle=\frac{1}{\sqrt{2}^{2n}}\langle\mathbf{x}^{\sigma}\mathbf{x}^{\tau}|H^{\otimes