Quantum supremacy in driven quantum many-body systems

A crucial milestone in the field of quantum simulation and computation is to demonstrate that a quantum device can compute certain tasks that are impossible to reproduce by a classical computer with any reasonable resources. Such a demonstration is referred to as quantum supremacy. One of the most important questions is to identify setups that exhibit quantum supremacy and can be implemented with current quantum technology. The two standard candidates are boson sampling and random quantum circuits. Here, we show that quantum supremacy can be obtained in generic periodically-driven quantum many-body systems. Our analysis is based on the eigenstate thermalization hypothesis and strongly-held conjectures in complexity theory. To illustrate our work, We give examples of simple disordered Ising chains driven by global magnetic fields and Bose-Hubbard chains with modulated hoppings. Our proposal opens the way for a large class of quantum platforms to demonstrate and benchmark quantum supremacy.

Introduction-Quantum computational supremacy is the ability of quantum devices to efficiently perform certain tasks that cannot be efficiently done on a classical computer [1,2]. Early proposals for realizing quantum supremacy include boson sampling [3][4][5] and random quantum circuits [6][7][8]. In both cases, the computational hardness stems from the inability of a classical computer to efficiently sample the output probabilities of a complex quantum evolution. Experimental efforts towards achieving quantum supremacy include optical networks for boson sampling [9][10][11][12][13] and superconducting circuits for random circuits [14]. Signatures of quantum supremacy have been observed recently with 53 superconducting qubits [15].
Analog quantum simulators are controllable quantum platforms specifically built to implement complex quantum many body models [16][17][18][19]. In these experiments, complex quantum dynamics have been implemented which cannot be reproduced with existing classical numerics and have shed light on important questions in quantum many-body physics [20]. However, rigorous proof of quantum supremacy involving complexity theory in those systems are yet to be shown, with the few exceptions of the 2D quantum Ising [21,22] and the 2D cluster-state models [23].
In this work, we provide evidence that when generic isolated periodically-driven quantum many-body systems thermalize, in the sense that any observables can be obtained from the microcanonical ensemble, sampling from their output distribution cannot be efficiently performed on a classical computer. These constitute a large class of quantum simulators that are currently available [14,[24][25][26][27][28][29]. Our analysis is based on the absence of collapse of the polynomial hierarchy and two plausible assumptions: the worst-to average-case hardness of the sampling task and the experimental realisability of random * cqtjt@nus.edu.sg † dimitris.angelakis@qubit.org unitaries as predicted by the Floquet Eigenstate Thermalization Hypothesis (ETH). We support our findings by examining specific examples of disordered quantum Ising chain driven by a global magnetic field and the onedimensional Bose-Hubbard (BH) model with modulated hoppings. These models have been widely implemented experimentally [14,[24][25][26][27][28][29], making our work of broad interest to the experimental community. General framework-Let us consider a generic periodically-driven quantum many-body system whose Hamiltonian is described byĤ(t) =Ĥ 0 +f (t)V . HereĤ 0 is the undriven Hamiltonian,V is the driving Hamiltonian such that Ĥ 0 ,V = 0, and f (t) is periodic with period T . We require that the time-averaged Hamiltonian H ave = 1 T T 0Ĥ (t)dt describes an interacting many-body system [30].
Let Z = {|z = ⊗ L i |z i } be a complete basis of manybody Fock states, where z i = {0, 1, 2, .., D i − 1} denotes the basis state of a local quantum system of dimension D i and where i ∈ [1, L]. In what follows, we assume without loss of generality that D i = D for all i, resulting in an Hilbert space of dimension N = D L . The state after M driving periods is |ψ M =Û M F |z 0 , whereÛ F = T exp −i T 0Ĥ (t)dt ≡ exp −iĤ F T andT is the time-ordering operator. We assume that the initial state |z 0 is a product state. The effective time-independent Floquet HamiltonianĤ F fully describes the dynamics probed at stroboscopic times t = nT . The probability of measuring the Fock state |z is then p M (z) = | z|ψ M | 2 with where the sum is performed over M − 1 complete sets of basis states. More precisely, the set of basis states {|z m } is associated with the quantum evolution after m driving cycles with z 0 (z M = z) being the initial (readout) configuration. The expression in Eq. (1) can be viewed as the Feynman's path integral where each trajectory is defined by a set of configurations {z 0 , z 1 , ..., z M }. The ETH states that generic isolated many-body quantum systems thermalize by their own dynamics after a long enough time, regardless of their initial state. In that case, any generic observable is expected to evolve toward the canonical ensemble with a finite temperature [31]. For driven quantum many-body systems, it has been shown that not only thermalization still occurs, but that for low-frequency driving, the associated temperature becomes infinite [32]. In this limit, the Floquet operatorÛ F shares the statistical properties of the Circular Orthogonal Ensemble (COE). This is an ensemble of matrices whose elements are independent normal complex random variables subjected to the orthogonality and the unitary constraints. This emergent randomness is the particular ingredient responsible for the hardness in calculating the output probability of Eq. (1), as there are exponentially many random Feynman trajectories that are equally important. We emphasize that the external periodic drive is crucial to reach the required level of randomness [33,34]. A more detailed analysis ofĤ F shows that the presence of low-frequency driving allows to generate effective infiniterange multi-body interactions [32,35]. Therefore lifting the constraints imposed by the limited local few-body interactions generally encountered in physical systems [36].
Quantum supremacy-To understand the computational task, let us first define some essential terms used in the complexity theory, namely approximating, sampling, multiplicative error and additive error. Let us imagine an analog quantum device built to mimic the quantum dynamics that would lead to p M (z) = | z|ψ M | 2 . In practice, such device will encode an output probability q(z) that differs from p M (z) due to noise, decoherence and imperfect controls. Both probabilities are said to be multiplicatively close if where α ≥ 0. The task of approximating p M (z) up to multiplicative error is to calculate q(z) that satisfies the above equation for a given z. However, such degree of precision is difficult to achieve experimentally as the allowed error is proportional to p M (z) which can be much smaller than unity. A more feasible task is to approximate p M (z) up to additive error, defined as with β > 0. Note that the additive error involves summing over all possible output strings z ∈ Z, while the multiplicative condition applies to each z individually. The task of approximating p M (z) even with additive error is still unrealistic as it requires a number of measurements that grows exponentially with the size of the system. What a quantum device can do is to sample strings from q(z). Hence, we define the task of sampling from p M (z) up to additive error as generating strings from q(z) while q(z) is additively close to p M (z). This task is our central focus to show quantum supremacy. We emphasize that it is different from "certifying quantum supremacy" [37] which consists of certifying if Eq. (3) holds.
To show that the above sampling task cannot be done efficiently by a classical computer, we follow the standard argument which proceeds as follows. Let us suppose that there is a classical machine C able to sample from p M (z) up to additive error and that the distribution of p M (z) anti-concentrates, i.e.
Here, the distribution is obtained from a set of unitary matrices {Û F } that are realizable experimentally. The Stockmeyer theorem states that, with the help of a NP oracle, that machine C can also approximate p M (z) up to multiplicative error for some outcomes z [39]. We emphasize that the sampling task is converted to the approximation task in this step. If the latter is #P-hard, then the existence of that machine C would imply the collapse of the polynomial hierarchy to the third level, which is strongly believed to be unlikely in computer science. Hence, assuming that the polynomial hierarchy does not collapse to the third level, we reach the conclusion that a classical machine C does not exist. The two fundamental conditions of the proof, that is the #P-hardness of approximating p M (z) up to multiplicative error and the anti-concentration of p M (z), can be more formally based on the two following theorems.  In theorem 1, we introduced the key notion of worstcase hardness of the entire set of COE matrices {Û COE }. This corresponds to the scenario where at least one instancep M (z), i.e. a single unitaryÛ ∈ {Û COE } and a single configuration z ∈ Z, is hard to approximate with multiplicative error. However, that one instance may be impractical to produce experimentally as the full set of COE matrices {Û COE } (p M (z)) might not coincide with the experimentally accessible set {Û F } (p M (z)). That is due to the fact that even though Floquet ETH strongly suggests thatÛ M F is an instance uniformly drawn from the {Û M COE }, not allÛ M COE matrices will be realizable by {Û M F }. More desirable is the average-case hardness where most instances are hard. Consequently, to ensure that the hard instance in Y can be found within {Û F } and thus prove quantum supremacy for realizable driven analog quantum systems, we further assume the following two conjectures.
Conjecture 1 (Average-case hardness) For any 1/2e fraction of Y, approximatingp M (z) up to multiplicative error with α = 1/4 + o(1) is as hard as the hardest instance. Here o(·) is the little-o notation. Informally, conjecture 1 assumes the worst-to-average case reduction in Y which is common in most quantum supremacy proposals [40]. Conjecture 2 connects the mathematically constructed COE to experimentally accessible driven analog quantum systems by stating that the ensemble {Û M F } is statistically equivalent to a set of instances uniformly drawn from {Û M COE }. This conjecture is supported by the observation that isolated systems evolving under genericÛ M F thermalize to infinite temperature resulting in fully random final quantum states [32][33][34] in experimentally relevant timescales [24,35]. Compared to existing quantum supremacy proposals, the reliance of the main theorem (see below) on conjecture 2 is not standard and may be seen as undesirable. But from our perspective, this conjecture makes a connection between computational complexity and the experimentally tested Floquet ETH that is applicable to a broad class of generic periodically-driven quantum systems as implemented in a variety of analog quantum simulators. Proving or disproving conjecture 2, either directly or indirectly by refutation of the main theorem while conjecture 1 holds true, is by itself of fundamental interest in physics.
The fraction used in conjecture 1 is chosen so that the approximate Haar random measure ensures that some hard instances in Y can be realized with {Û M F }. In combination of theorems 1 and 2, the two conjectures finally allow us to state the main theorem of this work.
Main Theorem Assuming conjectures 1 and 2, the ability to classically sample from p M (z) up to an additive error β = 1/8e for all unitary matrices in {Û F } implies the collapse of the polynomial hierarchy to the third level.
In what follows, we address in detail the proofs of theorems 1 and 2 while the detailed application of the standard Stockmeyer argument to prove the main theorem is provided in Appendix A #P hardness of simulating COE quantum dynamics-To prove theorem 1, we first notice that the COE is an ensemble of all orthogonal unitary matrices. This includes the well-known instantaneous quantum polynomial (IQP) circuitsÛ IQP =ĤẐĤ, whereĤ consists of Hadamard gates andẐ is an arbitrary (possibly nonlocal) diagonal gate on the computational basis, both acting on all qubits [6]. The IQP circuits constitute one of the early proposals of quantum supremacy. Multiplicative approximation of their output probabilities are known to be #P-hard in the worst case [41,Theorem 1.4]. SinceÛ M IQP =ĤẐ MĤ still adopt the general form of the IQP circuits, we conclude that there exists at least one instance in Y that is #P-hard for multiplicative approximation.
To see how the hardness could emerge for a typical instance in Y (conjecture 1), one can in principle map the path integral in Eq. (1) to the partition function of a classical Ising model with random complex fields. The latter is widely conjectured to be #P-hard on average for multiplicative approximation [21,42]. In this context, the key is to note that a COE unitary evolution can be written asÛ COE =Û T CUEÛ CUE , whereÛ CUE is a random matrix drawn from the Circular Unitary Ensemble (CUE), i.e. the ensemble of Haar-random matrices [43]. Furthermore,Û CUE can be decomposed into a set of universal quantum gates which can be mapped onto a complex Ising model. This mapping procedure has already been described in ref. [7] to support the conjecture of the worst-to-average case in the context of random quantum circuits. A detailed and intuitive description of this protocol is presented in Appendix B.
Anti-concentration of COE dynamics.-To prove the second and necessary ingredient of the proof, i.e. theorem 2, we write where d (z) = z|E E |z 0 , φ M, = M E T mod 2π, |E is an eigenstate ofĤ F with eigenenergy E . For COE operators, d (z) are real [43] and their distribution, denoted as Pr(d), is given by the Bessel function of the second kind (see Fig. 1(a) and Appendix C for a detailed derivation). Consequently, the values of d (z) for different and z do not concentrate on a particular value. Now let us consider the statistics of the phases {φ M, }.
We define the level spacing as For a single driving cycle M = 1, the phases {φ 1, } for COE are known to exhibit phase repulsion, i.e. the phases are correlated [32]. The COE distribution Pr COE (r ) is depicted in Fig. 1(b), where Pr COE (0) = 0 explicitly indicates the phase repulsion. For multiple driving cycles M 2π/E T , the correlations are erased due to energy folding, i.e. the effect of the modulo 2π. This results in the Poisson (POI) distribution of the level spacing, Pr POI (r ) = 2/(1 + r 2 ), with the peak at r = 0, see Fig. 1 The Bessel function distribution of d (z) and the POI distribution of φ M, ensure that the output distribution Pr which suggests that the system explores uniformly (approximately Haar-random) the Hilbert space [14,21]. This satisfies the anti-concentration condition since [7]. To see the emergence of the Porter-Thomas distribution, we write Due to the Poisson distribution in the long time limit, the phases {φ M, } can be thought of as independent variables randomly and uniformly distributed in the range [0, 2π). Using the product distribution formula and the central limit theorem, one can show that the distributions of a z and b z are normal distributions with zero mean and variance 1/2N . Sincep M (z) = a 2 z + b 2 z , the Porter-Thomas distribution ofp M (z) can be derived using the fact that the square sum of two Gaussian variables follows the χ-squared distribution with second degree of freedom [44]. A detailed derivation is presented in Appendix C.
Example of driven many-body systems.-We give two specific examples of driven systems that display statistical properties consistent with the COE and hence partially support conjecture 2. For both cases, the modulation is f (t) = 1 2 (1 − cos(ωt)), where ω = 2π/T and initial states are randomized product states.
(i) 1D Ising chain: We consider an Ising chain described by the HamiltonianĤ ISING W is the disorder strength,Ẑ l is the Pauli spin operator acting on site l, and J is the interaction strength. The drive is a global magnetic fieldV ISING = F L−1 l=0X l , where F is the driving amplitude. Similar models have been implemented in various quantum platforms, including trapped ions [27] and superconducting circuits [28].
(ii) 1D Bose-Hubbard model: We consider the

BH model described by the HamiltonianĤ
is a bosonic annihilation (creation) operator at site l, U is the onsite interaction, and µ l is the local disorder as defined above. The drive modulates the hopping amplitudeŝ ). Similar models have been implemented in superconducting circuits [14] and cold atoms [24,26,29].
The distribution of d (z) from both models are depicted in Fig. 1(a), showing an agreement with the Bessel function as predicted by COE. The level statistics at M = 1 and M = 25 are depicted in Fig. 1(b), showing an agreement with the COE and the POI distribution, respectively. The driving frequency and the disorder strength are tuned to ensure the observation of the thermalized phase and prevent many-body localization [32,45]. Fig. 2 shows the l 1 -norm distance between Pr(p) and the Porter-Thomas distribution at different m for the Ising and the BH models. It can be seen that, in all cases, the system reaches the Porter-Thomas distribution after multiple driving cycles. The l 1 -norm distance in the longtime limit is decaying towards zero as the size of the system increases. Therefore, the anti-concentration condition is satisfied. In absence of the drive, a similar analysis can be performed for the infinite-time unitary evolution corresponding to generic instances of the undriven thermalized phase in both models. In this case, d (z) does not follow the Bessel function of the second kind and the output distribution never reaches the Porter-Thomas distribution (see Appendix D for numerical simulation of the undriven Ising and Bose-Hubbard models). This is consequence of the energy conservation and the structure imposed by the local interactions, highlighting the key role played by the drive.
Conclusions and outlook-Analog quantum simulators realizing quantum many-body systems have generated quantum dynamics beyond the reach of existing classical numerical methods for some time. However, such dynamics has not been theoretically proven to be hard to compute by a classical computer. We have shown here that in the particular case of driven many-body systems, when they thermalize, sampling from their output distribution cannot be efficiently performed on a classical computer. Using complexity theory arguments, we provide strong analytical evidence of the computational hardness stemming from the COE statistics, and provide numerical results showing that COE dynamics can be obtained from driven quantum Ising and BH models for realistic parameters.
Our results greatly widen the possibilities to realise quantum supremacy with existing experimental platforms and provide the theoretical foundations needed to demonstrate quantum supremacy in analog quantum simulators. In the future, it would be interesting to extend our results to a broader class of quantum many-body systems such as those with gauge fields, frustrated spin systems, and undriven systems. For example, in Ref. [20], cold atoms in optical lattices have been used to compute the undriven quantum many-body localization transition in two dimensions, which has so far eluded state-of-theart classical numerical techniques [46].
In this section, we provide a detailed proof of the main theorem of the main text, which reads: Main Theorem Assuming conjecture 1 and 2, the ability to classically sample from p M (z) up to an additive error β = 1/(8e) for all unitary matrices in {Û F } implies the collapse of the polynomial hierarchy to the third level.
The proof relies on the theorems 1 and 2 and conjectures 1 and 2 presented in the main text.   (1), where o(·) is little-o notation [49], is as hard as the hardest instance. Let us begin by considering a classical probabilistic computer with an NP oracle, also called a BPP NP machine. This is a theoretical object that can solve problems in the third level of the polynomial hierarchy. The Stockmeyer theorem states that a BPP NP machine with an access to a classical sampler C, as defined in the main text, can efficiently output an approximationq(z) of q(z) such that |q(z) −q(z)| ≤ q(z) poly(L) .
We emphasise that the BPP NP machine grants us the ability to perform the approximating task, in contrast to the machine C that can only sample strings from a given distribution. To see how the BPP NP machine can output a multiplicative approximation of p M (z) for most of z ∈ Z, let us consider .
The first and the third lines are obtained using the triangular inequality. To get multiplicative approximation of p M (z) usingq(z), we need the term |p M (z) − q(z)| to be small. Given the additive error defined in Eq. (3) in the main text, this is indeed the case for a large portion of {z} ∈ Z. Since the left hand side of Eq. (3) in the main text involves summing over an exponentially large number of terms but the total error is bounded by a constant β, most of the terms in the sum must be exponentially small. This statement can be made precise using Markov's inequality.

Fact 1 (Markov's inequality)
If X is a non-negative random variable and a > 0, then the probability that X is at least a is where E(X) is the expectation value of X.
Here, the distribution and the expectation value are computed over z ∈ Z. Note that E z (|p M (z) − q(z)|) ≤ β/N is given by the additive error defined in Eq. (3) in the main text. By setting a = β/N ζ for some small ζ > 0, we get By substituting |p M (z) − q(z)| from Eq. (A2), we get Theorem 2 in the main text (the anti-concentration condition) together with conjecture 2 imply that {p M (z)} follows the Porter-Thomas distribution, specially that 1/N < p M (z) for at least 1/e fraction of the unitary matrices in {Û F }. Hence, we can rewrite Eq. (A7) as Here, the distribution is over all z ∈ Z and all unitary matrices in {Û F }. To understand the right hand side of the equation, let P ∩ Q be the intersection between the set P of probabilities that anticoncentrate and the set Q of probabilities that satisfy the Markov's inequality. Since Pr(P ∩ Q) = Pr(P ) + Pr(Q) − Pr(P ∪ Q) ≥ Pr(P ) + Pr(Q) − 1, Pr(P ) = 1/e and Pr(Q) = 1 − ζ , it follows that Pr(P ∩ Q) is no less than 1/e + 1 − ζ − 1 = 1/e − ζ.
Following [7,42], we further set β = 1/(8e) and ζ = 1/(2e), so that giving an approximation up to multiplicative error 1/4 + o(1) for at least 1/(2e) instances of the set of experimentally realizable unitary matrices {Û F }. If according to the conjecture 1 and conjecture 2 in the main text, multiplicatively estimating 1/(2e) fraction of the output probabilities from {Û F } is #P-hard, then the Polynomial Hierarchy collapses. This concludes the proof of the main theorem in the main text.
Appendix B: Mapping of approximating output distribution of COE dynamics onto estimating partition function of complex Ising models In this section, we provide evidence to support the conjecture 1 in the main text, showing how hardness instances could appear on average. To do this, we map the task of approximating an output distributions of COE dynamics onto calculating the partition function of a classical Ising model which is widely believed to be #P-hard on average for multiplicative approximation [21,42]. The section is divided into two parts. In the first part, we explain the overall concept and physical intuition of this procedure. In the second part, mathematical details are provided.

Physical perspective of the mapping procedure
The mapping protocol consists of two intermediate procedures. First, we map the COE unitary evolution on universal random quantum circuits and, second, we derive a complex Ising model from those circuits following Ref. [7].
Let us begin by expressing an unitary evolution of COE asÛ COE =Û T CUEÛ CUE whereÛ CUE is a random unitary drawn from the Circular Unitary Ensemble (CUE) i.e. Haar ensemble [43]. We then further decomposeÛ CUE into a set of universal quantum gates [7]. Following Ref. [7], we choose random quantum circuits consisting of n + 1 layers of gates and log 2 N qubits, as shown in Fig. 3(a). The first layer consists of Hadamard gates applied to all qubits. The following layers consist of randomly chosen single-qubit gates from the set { √ X, √ Y , T } and two-qubit controlled-Z (CZ) gates. Here, √ X ( √ Y ) represents a π/2 rotation around the X (Y ) axis of the Bloch sphere andT is a non-Clifford gate representing a diagonal matrix {1, e iπ/4 }. Such circuits have been shown to be approximately t-design [50] for an arbitrary large t when n → ∞, which implies the CUE evolution [51]. The operatorÛ T CUE can be implemented by reversing the order of the gates inÛ CUE and replacing √ Y with √ Y T . We emphasize that decomposing the COE evolution into the random circuits is only done theoretically with an aim to show the average case hardness. In the real experiments, this COE dynamics is realized by the driven many-body systems. The mathematical procedure for the mapping from random quantum circuits to classical complex Ising models is discussed in details in the next part. Specifically, p M (z) from the circuit (Û T CUEÛ CUE ) M , as depicted in Fig. 3(a), can be calculated from the partition function, Here, A(s) is the degeneracy number associated with a classical spin configuration s in the lattice S, s i = ±1, h i represents a on-site field on site i and J ij represents the coupling between the classical spins on site i and j. Since the output probability can also be interpreted as the path integral in Eq. (B1) in the main text, the intuition behind the mapping is that the sum over all possible paths is translated into the sum over all possible classical spin configurations, where the phase accumulated in each path is given by the energy of the complex Ising lattice S. To gain intuitive understanding of this standard mapping, we provide a diagrammatic approach to visualize the lattice S and extract the field parameters {h i }, {J ij }. To begin with, we use the random circuit in Fig. 3(b) as a demonstration. The mathematical descriptions behind each steps are discussed in the next part.
• STEP I -For each qubit, draw a circle between every consecutive non-diagonal gates, see Fig. 3(c). Each circle or 'node' represents one classical spin. • STEP II -For each qubit, draw a horizontal line between every consecutive nodes i,j, see Fig. 3(d). These lines or 'edges' represent interaction J ij between two neighboring spins in the same row. In addition, draw a line between every two nodes that are connected by CZ gates. These lines represent the interaction J ij between spins in different rows.
• STEP III -Labeling each nodes and edges with the corresponding gates, see Fig. 3(e).
• STEP IV -Use the lookup table in Fig. 3(f) to specify h i and J ij introduced by each gate. For example, the √ Y gate that acts between nodes i and j adds −1 to J ij , −1 to h i and +1 to h j . We use the convention that the leftmost index represents the leftmost node. Also, the two T-gates that are enclosed by the node i will add 0.5 + 0.5 = +1 to the local field h i .
• STEP V -Finally, spins at the leftmost side of the lattice are fixed at +1, corresponding to the initial state |0 .
Similarly, spins at the rightmost side of the lattice are fixed according to the readout state |z .
Following the above recipe, we provide the exact form of the parameters in the Ising model for the COE dynamics in the next part, showing that the field parameters {h i } and {J ij } are quasi-random numbers with no apparent structure. Specifically, neither the phase π i h i s i /4 nor the phase π i,j J ij s i s j /4 is restricted to the values 0, π/2, π, 3π/2 (mod 2π) for each spin configurations s. Without such stringent restrictions, approximating the partition function up to multiplicative error is known to be #P-hard in the worst case [41,Theorem 1.9]. This motivates a widely used conjecture in quantum supremacy proposals that such task is also hard on average [21,42]. We emphasize here the major differences between random quantum circuits as proposed in Ref. [7] and our systems. Firstly, our systems are analog with no physical quantum gates involved. The decomposition to quantum gates is only done mathematically. Secondly, our system has discrete time-reversal symmetry, while such symmetry is absent in random quantum circuits. Consequently, the COE in our system is achieved from the Floquet operatorÛ F , while the CUE in random quantum circuits are achieved from the entire unitary evolution. In addition,Û M F in our system does not have the t-design property due to the COE [52, pp.117-119]. However, as shown above, the hardness arguments for the random quantum circuits can be naturally applied to our case.

Mathematical details of the mapping procedure
In this section, we prove Eq. (B1) by providing justifications of the diagrammatic recipes to map the the evolution U CUE on a Ising spin model with complex fields. Again, the quantum gates of interest consist of both diagonal gates {T, CZ} and non-diagonal gates For simplicity, we start with one-and two-qubit examples before generalizing to the COE dynamics. The mathematical procedure here is adapted from Ref. [7].

a. One-qubit example
Let us consider a one-qubit circuit and N + 1 gates randomly chosen from the set gate is fixed to be a Hadamard gate. The output probability is p(z) = | z|Û |0 | 2 , whereÛ = N n=0Û (n) is the total unitary matrix,Û (n) is the n th gate and z ∈ {0, 1} is the readout bit. Below, we outline the mathematical steps underlying the diagrammatic approach followed by detailed explanations for each step: In the second line, we insert an identityÎ n = zn∈{0,1} |z n z n | betweenÛ (n+1) andÛ (n) for every n ∈ {0, .., N − 1}. As a result, this line can be interpreted as the Feynman's path integral where each individual path or 'world-line' is characterized by a sequence of basis variables z = (z −1 , z 0 , ..., z N ). The initial and the end points for every path are |z −1 = |0 and |z N = |z , respectively. In the third line, we decompose z n |Û (n) |z n−1 into the amplitude A(z n , z n−1 ) and phase Φ(z n , z n−1 ). In the fourth line, we introduce A(z) = N n=0 A(z n , z n−1 ). The equation now takes the form of the partition of a classical Ising model with complex energies. Here, z can be interpreted as a classical spin configuration, A(z) as the degeneracy number and i π 4 Φ(z n , z n−1 ) as a complex energy associated with spin-spin interaction.
Further simplifications are possible by noticing that, the diagonal gates in the circuits allow the reduction of the number of classical spins. Specifically, if a T gate is applied to |z n−1 , it follows that z n = z n−1 . Hence, the variables z n−1 and z n can be represented by a single classical spin state. The two variables z n−1 , z n become independent only when a non-diagonal gate is applied. Therefore, we can group all variables {z n } between two non-diagonal gates as one classical spin. This procedure leads to the directives presented as the the STEP I of the procedure in the previous section. Formally, for N spin + 1 non-diagonal gates in the circuit (including the first Hadamard gate) z can be characterized by a classical spin configuration s = (s −1 , s 0 , ..., s k , ..., s Nspin ) where s k = 1 − 2z k ∈ {±1} is a spin representing the basis variable immediately after the k th non-diagonal gate, i.e.
Lastly, we need to specify A(s) and Φ(s k , s k−1 ) in term of the local fields h k−1 , h k , the interaction J k−1,k , and spin configurations s k−1 , s k . This is done by first considering the gates in their matrix form, i.e.
Notice that all non-diagonal gates contribute to the same amplitude A(s k , s k−1 ) = 1/ √ 2, leading to A(s) = 2 −(Nspin+1)/2 . Hence, we can extract the contribution of each gate to Φ(s k , s k−1 ) as The under-script indicates which gate is contributing to the phase. The corresponding h i , h j and J ij are depicted in the lookup table in Fig. 3(f), where i = k − 1 and j = k. The global phase that does not depend on s is ignored as it does not contribute to p(z).

b. Two-qubit example
Now we consider a two-qubit random circuits to demonstrate the action of the CZ gates. We introduce a new index l ∈ {1, 2} to label each qubit, which is placed on a given horizontal line (row). Since the CZ gate is diagonal, its presence does not alter the number of spins in each row. However, the gate introduces interaction between spins in different rows. This can be seen from its explicit form, i.e.
where s 1,k (s 2,k ) is the state of the k th (k th ) spin at the first (second) row. It follows that The corresponding h i , h j , and J ij are depicted in Fig. 3(f) where i = (1, k) and j = (2, k ). We have now derived all necessary ingredients to map a random quantum circuit to a classical Ising model.

c. Full COE dynamics
Since the COE dynamics can be expressed in terms of a quasi-random quantum circuit, we can straightforwardly apply the above procedure to find the corresponding Ising model. The complexity here solely arises from the number of indices required to specify the positions of all the gates in the circuit. To deal with this, we introduce the following indices -an index l ∈ {1, ..., L} to indicate which qubit / row.
-an index µ ∈ {A, B} to indicate which part of the period. A and B refer to theÛ CUE part and theÛ T CUE part, respectively -an index k ∈ {0, 1, ..., N spin (l)} to indicate the spin position for a given m and µ. Here, N spin (l) is the total number of spins at the l th row. Note that due to the symmetric structure ofÛ CUE andÛ T CUE , we run the index k backward for the transpose part, i.e. k = 0 refers to the last layer.
-an index ν l,k so that ν l,k = 1 if the k th non-diagonal gate acting on the qubit l is √ X otherwise ν l,k = 0.
With these indices, the partition function of the circuit, as shown in Fig. 3(a), can be written as Here G is the total number of non-diagonal gates in the circuit. ζ  In this section, we provide additional mathematical details involved in the proof of theorem 2. More precisely, we show that the distribution of the output probability of COE dynamics, Pr(p), follows the Porter-Thomas distribution Pr PT (p) = N e −N p . First, let us consider the output probability p M (z) = | z|ψ M | 2 with To prove lemma 1, we first write d (z) = c z, c 0, , where c z, = z|E and c 0, = 0|E . For the COE dynamics, the coefficients c z, and c 0, are real numbers whose distribution is [43] As discussed in the main text, the phase φ M, becomes random as M 2π/E T . The random sign (±1) from c z, can therefore be absorbed into the phase without changing its statistics. The distribution of d (z) can be obtained using the product distribution formula where K 0 is the modified Bessel function of the second kind.
To prove lemma 2, we first note that the distribution of cos φ m, and sin φ m, with φ M, being uniformly distributed in the range [0, 2π) are Pr(cos φ) = 1 Pr(sin φ) = 1 We then calculate the distribution of κ ≡ d (z) cos φ M, using the product distribution formula, i.e. The mean and the variance of κ can be calculated as Since a z is a sum of independent and identically distributed random variables, i.e. a z = N −1 =1 κ , we can apply the central limit theorem for large N . Hence, the distribution of a z is normal with the mean zero and variance Var(a) = N · Var(κ) = 1/2N . The same applies for the distribution of b z . Theorem 2 can be proven using the fact that the sum of the square of Gaussian variables follows the χ-squared distribution with second degree of freedom Pr χ 2 ,k=2 (p) ∼ exp{−p/2σ 2 } [44]. By specifying the variance obtained in Lemma 2 and normalization, the distribution of p M (z) = a 2 z + b 2 z over ∀z ∈ {0, 1} L is the Porter-Thomas distribution. Since the Porter Thomas distribution anti-concentrates i.e. Pr PT p > 1 N = ∞ N p=1 d(N p)e −N p = 1/e , we complete the proof of the theorem 2.
Appendix D: Undriven thermalized many-body systems In this section, we analyze the long-time unitary evolution for undriven systems in the thermalized phase. The results presented here highlight the key role played by the drive in generating the randomness required for the above quantum supremacy proof. In particular, we show that for typical undriven physical systems with local constraints (e.g. finite-range interactions) and conserved energy, the output distribution never coincides with the PT distribution.
We emphasize that this is a consequence of the inability of random matrix theory to accurately describe the full spectral range of undriven thermalized many-body systems. Indeed, it has been shown that for undriven many-body systems which thermalizes (to a finite temperature), the statistics of the Hamiltonian resembles the statistics of the Gaussian orthogonal ensemble (GOE) [31]. However, it is implicit that an accurate match only applies over a small energy window (usually far from the edges of the spectrum). If one zooms in this small energy window, the Hamiltonian looks random, but if one consider the full spectrum, the local structure of the Hamiltonian appears and the random matrix theory fails at capturing it.
To see this, we numerically simulate the undriven Ising Hamiltonian,Ĥ 0 = L−1 l=0 µ lẐl + J L−2 l=0Ẑ lẐl+1 + F 2 L−1 l=0X l , where µ l ∈ {0, W } is a local disorder, W is the disorder strength, F is the static global magnetic field along x and J is the interaction strength. This Hamiltonian is in fact the average Hamiltonian of the driven Ising Hamiltonian used in the main text. In comparison, we also simulate the quantum evolution under an ensemble {Ĥ COE } of synthetic Hamiltonians that are uniformly drawn from the GOE (i.e. without any local constraints). Fig.4 (a) shows the level-spacing statistics of {Ĥ 0 } (obtained over 500 disorder realizations), {Ĥ COE } (obtained over 500 random instances) and their corresponding long-time unitary operatorsÛ = lim t→∞ e −itĤ . We see that the level statistic of the physical Hamiltonian (and its long-time evolution) is indistinguishable from the GOE. However, the discrepancy between the physical and synthetic (GOE) realizations becomes apparent when looking at the eigenstate statistics as shown in Fig.4 (b). While the distribution of d (z) [see Eq. (5) of the main text] from the GOE is in a good agreement with the Bessel function of the second kind, the physical system fails to meet the theoretical prediction. This is in contrast to the driven case as presented in the main text. More importantly in the context of this work, a key difference between the physical Hamiltonian and the random matrix theory prediction can be seen by comparing the distribution of the output states after some time evolution. In Fig.4 (c), we show that the Porter-Thomas distribution is never achieved with the physical systems while it is for the synthetic realizations as well as for the driven case studied in the main text. These results underline the gap between physical Hamiltonians and true random matrices and more importantly, they highlights the important role of the drive in bridging that gap.