# Quantum supremacy and quantum phase transitions

Supanut Thanasilp<EMAIL_ADDRESS>Centre for Quantum
Technologies, National University of Singapore, 3 Science Drive 2, Singapore
117543  Jirawat Tangpanitanon  Centre for Quantum Technologies, National
University of Singapore, 3 Science Drive 2, Singapore 117543  Marc-Antoine
Lemonde  Centre for Quantum Technologies, National University of Singapore, 3
Science Drive 2, Singapore 117543  Ninnat Dangniam  Department of Physics and
Center for Field Theory and Particle Physics, Fudan University, Shanghai
200433, China  Dimitris G. Angelakis<EMAIL_ADDRESS>Centre for
Quantum Technologies, National University of Singapore, 3 Science Drive 2,
Singapore 117543  School of Electrical and Computer Engineering, Technical
University of Crete, Chania, Greece 73100

###### Abstract

Demonstrating the ability of existing quantum platforms to perform certain
computational tasks intractable to classical computers represents a
cornerstone in quantum computing. Despite the growing number of such proposed
“quantum supreme” tasks, it remains an important challenge to identify their
direct applications. In this work, we describe how the approach proposed in
Ref. [arXiv:2002.11946] for demonstrating quantum supremacy in generic driven
analog many-body systems, such as those found in cold atom and ion setups, can
be extended to explore dynamical quantum phase transitions. We show how key
quantum supremacy signatures, such as the distance between the output
distribution and the expected Porter Thomas distribution at the supremacy
regime, can be used as effective order parameters. We apply this approach to a
periodically driven disordered 1D Ising model and show that we can accurately
capture the transition between the driven thermalized and many-body localized
phases. This approach also captures the transition towards the Floquet
prethermalized regime for high-frequency driving. Revisiting quantum phases of
matter under the light of the recent discussions about quantum supremacy draws
a link between complexity theory and analog many-body systems.

##  I Introduction

Recent experimental developments [1, 2, 3, 4, 5] have made the first steps
toward demonstrating the ability of noisy intermediate-scale quantum (NISQ)
devices [6] to perform certain computational tasks intractable to classical
computers [7, 8, 9, 10, 11, 12, 13, 14, 5, 15, 16, 17]. Early proposals to
demonstrate this so-called quantum supremacy include sampling from random
superconducting circuits [10, 11, 18] and boson sampling on photonic platforms
[13, 14, 5]. Despite being an important milestone for quantum computation in
the NISQ era, extending the proposed computational tasks to address real-world
problems remains an open challenge. This is particularly true in the context
of sampling from quantum random circuits where any practical applications
remains elusive. For example, the randomness required in such circuits makes
them in principle unsuitable for medium- to large-scale variational quantum
algorithms [19], thus greatly limiting their applications in near-term quantum
chemistry [20, 21], material sciences [22] and quantum machine learning [23].

More recently, the computational task of sampling from generic closed
periodically-driven quantum many-body systems has also been proposed to
demonstrate quantum supremacy [24]. This opens the road for a plethora of
analog based quantum simulators on periodically driven cold atoms and ions to
be used for such tasks [25, 26]. The computational complexity of sampling from
such dynamics is intimately related to the ability of the closed system to
thermalize under the combined influence of interactions and the external
drive. Once thermalised, the associated temperature is effectively infinite
resulting in a quantum evolution capable of exploring its entire Hilbert space
[27, 28, 29, 30]. This is in contrast with the driven many-body localized
(MBL) phase where the presence of disorder prevents the system to thermalize
and limits the quantum evolution to a restricted portion of the Hilbert space
[30, 29, 31, 32, 33, 34]. Efficient theoretical representations exploiting
this fact suggests that the sampling from MBL dynamics is tractable to
classical computers [35, 36, 37]. Therefore, the experiment we proposed in
Ref. [24] to achieve quantum supremacy can be interpreted as a protocol to
probe how much of the Hilbert space has been explored under specific driven
quantum dynamics. This ability of tuning in and out of the chaotic evolution
underlying the thermal phase by controlling the level of disorder has already
been exploited as a potential application in the context of quantum machine
learning [38].

In this work, we describe how the quantum supremacy experiment proposed in
Ref. [24] can be extended to probe the phase diagram of periodically-driven
disordered many-body systems. More specifically, we propose utilising key
quantum supremacy signatures as order parameters, such as the Kullback-Leibler
divergence (KLD) [39] between the output probability distribution and the
Porter-Thomas (PT) distribution [40]. The one to one matching in the behaviour
of these two distributions is a strong indication that the system is capable
of exploring its entire Hilbert space [10] and thus intimately related to the
computational complexity of the sampling task 111This is similar to the growth
of the bipartite entanglement entropy in closed systems which also captures
the increase in the capacity of the system to explore its Hilbert space [45].
While a large entanglement entropy implies a classically intractable matrix
product state (MPS) representation [66], a vanishing KLD from the PT
distribution has been shown to be a key condition to prove that any classical
approach is intractable [10, 24].. Therefore, such order parameter offers
additional insights into the dynamics compared to the measure of the
population imbalance commonly used in experiments [34, 42, 43, 44]. In
addition, the measurement of the KLD from the PT distribution, performed only
in the computational basis, should also be experimentally more accessible
compared to probing the level-statistics [45] which is already challenging for
undriven systems [46] with no known efficient extension to the driven cases.
We show that probing the KLD from the PT distribution allows to distinguish
the driven thermalized phase from the driven MBL phase as well as from the
Floquet prethermalized regime [28, 47]. The latter corresponds to the scenario
where the drive frequency exceeds all relevant energy scales of the system,
leading to long-lived prethermalization dynamics. Although we focus on the KLD
between the output and the PT distributions as it is a straightforward quantum
supremacy signature, we emphasize that our work could be extended to any
suitable metrics associated with the computational complexity of sampling from
the different phases. For example, we discuss in Appendix B how the
entanglement entropy over multiple sub-systems can also be used as an
alternative supremacy signature for probing the phase diagram.

As a specific model, we focus on an isolated Ising spin chain with nearest-
neighbour interactions, disordered on-site energies and a periodically time-
varying global magnetic field. We analyse the dynamics using the Floquet
formalism which allows to highlight the role played by the effective long-
range multi-body interactions generated by the drive in the context of quantum
supremacy. Our work thus proposes a direct application of a quantum supremacy
experiment and present additional physical intuitions about the origin of the
computational complexity associated with sampling from a driven thermalized
quantum system.

Figure 1: Schematic of a periodically driven quantum Ising chain with nearest-
neighbor interactions and disordered onsite energies.

##  II Quantum phases in driven disordered Ising spin chains

###  II.1 Model

The particular model we consider is sketched in Fig. 1 and consists of a
periodically-driven Ising chain with disordered onsite energies. The
corresponding Hamiltonian reads $\hat{H}(t)=\hat{H}_{0}+\hat{H}_{d}(t)$, where

$\displaystyle\hat{H}_{0}=\sum_{i=1}^{L}h_{i}\hat{\sigma}^{z}_{i}+B_{0}\sum_{i=1}^{L}\hat{\sigma}_{i}^{x}+J\sum_{i=1}^{L-1}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{i+1},$
(1) $\displaystyle\hat{H}_{d}(t)=\delta B\cos{(\omega
t)}\sum_{i=1}^{L}\hat{\sigma}^{x}_{i}.$ (2)

Here $\hat{\sigma}_{i}^{x}$ ($\hat{\sigma}_{i}^{z}$) is the $x$ ($z$) Pauli
matrix acting on the $i^{\rm th}$ spin of an open chain of $L$ total sites,
$J$ represents the nearest-neighbour interaction strength and
$B(t)=B_{0}+\delta B\cos(\omega t)$ represents a periodic magnetic field along
the $x$-axis of frequency $\omega\equiv 2\pi/T$. We consider disordered static
magnetic fields along the $z$-axis with strengths
$h_{i}\in\left[-W/2,W/2\right]$ being drawn from a uniform distribution
ranging from $-W/2$ to $W/2$. Importantly, for finite disorder $W\neq 0$ the
Hamiltonian at different times $t_{1}$ and $t_{2}\neq t_{1}+nT$ with
$n\in\mathbb{Z}$ do not commute, i.e. $[\hat{H}(t_{1}),\hat{H}(t_{2})]\neq 0$,
and the contribution from the drive averages to zero over a period, i.e.
$\frac{1}{T}\int_{0}^{T}\hat{H}_{d}(t)=0$.

The time evolution over a driving period is described by the unitary operator

$\hat{U}=\hat{\mathcal{T}}\exp\left[-i\int_{0}^{T}\hat{H}(t)dt\right]\equiv\exp\left[-i\hat{H}_{F}T\right],$
(3)

where $\hat{\mathcal{T}}$ is the time-ordering operator. The time-independent
Hamiltonian $\hat{H}_{F}$ is known as the Floquet Hamiltonian and fully
describes the dynamics at stroboscopic times $t_{n}=nT$. While a specific form
of the Ising model is studied in this work, the results obtained can be
generalized for other driven many-body systems such as Bose-Hubbard model or
different variations of the Ising spin model. All together, such driven
systems have been experimentally implemented in various quantum platforms
including superconducting circuits [2], cold atoms [48, 47, 25] and trapped
ions [26].

###  II.2 Phases and level-statistics

Depending on the disorder strength $W/J$ and driving frequency $\omega/J$, the
system can be in three distinct regimes [28, 34]. At low disorder strengths,
the system thermalizes under its own dynamics if given enough time and follows
the eigenstate thermalization hypothesis [28]. Under the assumption that the
system can be considered closed over the entire process, the constant energy
input provided by the drive leads to an effective temperature that is infinite
[27, 28, 29, 30]. The time-scale over which this thermalization process takes
place strongly depends on the drive frequency [28, 47]. For low-frequency
drives, the system can efficiently respond and rapidly reaches this infinite-
temperature limit corresponding to what is known as the driven thermalized
phase [27, 28, 29, 30]. In the case where the driving frequency exceeds all
relevant energy scales of the system however, the time required for
thermalizing can be greatly extended leading to a long-lived prethermalization
regime [28, 47]. Finally, in presence of large disorders, the system fails to
thermalize at any time and is said to be in the driven MBL phase [31, 33, 30,
32, 29]. The two phases and the prethermalized regime are depicted in Fig. 2
(a).

Figure 2: Level statistics of the driven disordered Ising spin chain. (a)
Average level spacing $\langle r\rangle_{\hat{U}}$ of the full unitary
$\hat{U}$ as a function of the disorder strength $W/J$ and the driving
frequency $\omega/J$ for $L=9$ and $B_{0}=-\delta B=1.25J$. The three limiting
statistics are indicated in white and the solid black line represents $\langle
r\rangle_{\hat{U}}=0.51$. (b) A cut along the frequency $\omega=8J$ [vertical
dashed line in (a)] shows $\langle r\rangle_{\hat{U}}$ as a function of $W/J$
and captures the phase transition from driven thermalized to the driven MBL.
(c) A cut along $W=4J$ [horizontal dashed line in (a)] captures the transition
as a function of $\omega/J$ from the driven thermalized phase to the
prethermalized regime. The red triangles represent the mean value $\langle
r\rangle_{\hat{U}_{0}}$ computed only considering the evolution under
$\hat{H}_{0}$ for a time $t=2\pi/\omega$. (d)-(f) Full distribution ${\rm
Pr}(r)$ for the driven thermalized phase ($W=4J,\omega=4.2J$) following the
COE, MBL phase ($W=30J,\omega=8J$) following POI and the prethermalization
regime ($W=4J,\omega=20.1J$) converging towards the GOE respectively. In (d),
the statistics of $\hat{U}$ and $\hat{U}_{0}$ drastically differ while they
converge in (f); contrasting the importance of the the multi-body long-range
interaction induced by the drive at different $\omega$. Each data point
results from $100$ disorder realizations.

One of the standard approach to distinguish the different phases is based on
the notion of level statistics of the unitary operator $\hat{U}$ [27, 31]. Let
$|\phi_{n}\rangle$ be an eigenstate of the Floquet Hamiltonian with eigenvalue
$\epsilon_{n}$, i.e.
$\hat{H}_{F}|\phi_{n}\rangle=\epsilon_{n}|\phi_{n}\rangle$, it follows that

$\hat{U}=\sum_{n}e^{i\theta_{n}}|\phi_{n}\rangle\langle\phi_{n}|,$ (4)

where $\theta_{n}=\epsilon_{n}T\text{ modulo 2}\pi$. The level statistics is
defined as the normalized distribution $\text{Pr}(r)$ of the level spacing

$r_{n}=\frac{\text{min}(\delta_{n},\delta_{n+1})}{\text{max}(\delta_{n},\delta_{n+1})},$
(5)

with $\delta_{n}=\theta_{n+1}-\theta_{n}$ and $\theta_{n+1}\geq\theta_{n}$.

In the driven thermalized phase, the level statistics of $\hat{U}$ coincides
with the circular orthogonal ensemble (COE) statistics, given by

$\displaystyle\text{Pr}_{\rm COE}(r)=$
$\displaystyle\frac{2}{3}\left[\frac{\sin\left(\frac{2\pi r}{r+1}\right)}{2\pi
r^{2}}+\frac{1}{(r+1)^{2}}+\frac{\sin\left(\frac{2\pi}{r+1}\right)}{2\pi}\right.$
$\displaystyle\left.-\frac{\cos\left(\frac{2\pi}{r+1}\right)}{r+1}-\frac{\cos\left(\frac{2\pi
r}{r+1}\right)}{r(r+1)}\right].$ (6)

This result means that, from a statistical point of view, a generic $\hat{U}$
associated with this phase is equivalent to a unitary matrix randomly sampled
from the complete ensemble of unitaries that conserves the time-reversal
symmetry [49]. The fact that Pr${}_{\rm COE}(0)=0$ indicates phase repulsion
and thus correlations between the different eigenstates. A COE evolution in
low-frequency driven systems is an indicator of thermalization to infinite
temperatures [27].

The scenario is completely different for systems in the driven MBL phase where
the level statistics is described by a Poisson distribution [33],

$\text{Pr}_{\rm POI}(r)=\frac{2}{(1+r)^{2}}.$ (7)

This time, $r=0$ corresponds to the peak of $\text{Pr}_{\rm POI}(r)$
indicating a high level of degeneracies in the MBL phase. This absence of
phase repulsion is intuitive as spatially distant localized eigenstates are
uncorrelated.

In the limit where the driving frequency exceeds all relevant energy scales of
the system, a slightly different interpretation of the level statistics is
required. This comes as all phases $\theta_{n}=\epsilon_{n}T<2\pi$ (for
$\omega>{\rm max}[\epsilon_{n}]$) become strictly linearly related to the
eigenenergies of $\hat{H}_{F}$ so that the statistics of $\hat{U}$ must
coincide with the statistics of $\hat{H}_{F}$. In addition, for
$\omega\rightarrow\infty$, $\hat{H}_{F}\rightarrow\hat{H}_{0}$ (see section
IV.3) and one should expect to recover the statistics of the undriven Ising
model as the frequency increases. Similar to the driven case, $\hat{H}_{0}$
has two distinct phases: the MBL phase for large disorder and the finite-
temperature thermalized phase for weak disorder. While the undriven MBL still
leads to the Poisson statistics, $\text{Pr}(r)$ of the undriven thermalized
phase follows the Gaussian orthogonal ensemble (GOE),

$\displaystyle\text{Pr}_{\rm
GOE}(r)=\frac{27}{4}\frac{r+r^{2}}{(1+r+r^{2})^{5/2}},$ (8)

which corresponds to the ensemble of matrices with independent normal random
variables as elements and subjected to the orthogonality constraint [49]. This
high frequency limit is known as the prethermalization regime and we note that
in the thermodynamics limit, the energy spectrum of the system becomes
infinite so that the driving frequency required to reach this regime also
becomes infinite.

In Fig. 2, we show the level statistics as a function of the drive frequency
$\omega/J$ and the disorder strength $W/J$ for an Ising chain of $L=9$. The
results are obtained by exactly diagonalizing $\hat{U}$ of Eq. (3) over 100
disorder realizations for each data point. In panels (a)-(c), we show the mean
value of level spacing $\langle r\rangle_{\hat{U}}$ as a metric to
characterize the different phases and compare it to the values obtained in the
COE ($\langle r\rangle_{\rm COE}\approx 0.527$), the POI ($\langle
r\rangle_{\rm POI}\approx 0.386$) and the GOE ($\langle r\rangle_{\rm
GOE}\approx 0.536$). We distinguish the three distinct regimes of the driven
dynamics and explicitly show the transition from driven thermalized to driven
MBL phases as we increase $W/J$ with fixed $\omega/J$ in panel (b) and from
driven thermalized phase to prethermalization regime as we increase $\omega/J$
with fixed $W/J$ in panel (c). In panels (d)-(f), we show Pr$(r)$ obtained
from $\hat{U}$ in these three regimes where it coincides with their respective
statistical ensembles. Finally, we also compute the level statistics of the
unitary operator produced by only considering the undriven Hamiltonian
$\hat{H}_{0}$, i.e. $\hat{U}_{0}=\exp{-i\hat{H}_{0}T}$, and compare it to the
exact results corresponding to the same frequency $\omega=2\pi/T$. As
expected, $\langle r\rangle_{\hat{U}}$ and $\langle r\rangle_{\hat{U}_{0}}$
converge as $\omega/J$ increases. We note that in absence of disorder,
$W\rightarrow 0$, the system evades thermalization for any driving
frequencies. This comes as the dynamics described by Eq. (1) with $W=0$ is
equivalent to a non-interacting fermionic system [50].

##  III Quantum supremacy in driven analog many-body systems

In this section, we discuss the computational complexity of sampling from
quantum states evolving under the dynamics associated with the different
regimes discussed above. We begin by briefly reviewing the formal evidences
presented in Ref. [24] which suggests that such a task is classically
intractable in the driven thermalized phase following the COE. We focus on the
convergence of the distribution of the output probabilities towards the PT
distribution as a key quantum supremacy signature which we contrast in the MBL
phase and the prethermalized regime.

###  III.1 Quantum supremacy in the driven thermalized phase

We define the exact probability of measuring the basis state $|{\bf
z}\rangle=\otimes_{i}^{L}|z_{i}\rangle$ after $m$ driving cycles as
$p_{m}({\bf z})\equiv|\langle{\bf z}|\psi_{m}\rangle|^{2}$, where
$\hat{\sigma}_{i}^{z}|z_{i}\rangle=z_{i}|z_{i}\rangle$ is the eigenstate of
the single-qubit Pauli $z$ matrix acting on site $i$ and where
$|\psi_{m}\rangle=(\hat{U})^{m}|\psi_{0}\rangle$. The initial state
$|\psi_{0}\rangle$ is assumed to be a product state. The computational task
considered here is the sampling of output bitstrings ${\bf z}$ from
$p_{m}({\bf z})$ up to (realistic) additive errors. More precisely, it
requires to produce a sampler $\mathcal{C}$ capable of providing bitstrings
from a distribution $q(\textbf{z})$ that is additively closed to the exact
distribution $p_{m}(\textbf{z})$, i.e.

$\sum_{\textbf{z}}|p_{m}(\textbf{z})-q(\textbf{z})|\leq\beta_{0},$ (9)

where $\beta_{0}$ is a small positive constant. The sampler $\mathcal{C}$ can
be either constructed by a quantum or a classical device. Quantum-
mechanically, this can be achieved by performing projective measurements on an
analog quantum simulator capable of reproducing the exact driven Hamiltonian.
The difference between $q(\textbf{z})$ and $p_{m}(\textbf{z})$ can be kept
small if the noise from the environment and imperfections in the controls
remains sufficiently low.

Producing this sampler $\mathcal{C}$ with a classical computer is however not
as natural and strong analytical evidence suggests that, given that $\hat{U}$
follows the COE statistics, it is classically intractable [24]. One of the key
features that ensures this computational complexity is the fact that the
associated quantum dynamics anti-concentrates [51]. Mathematically, this means
that the distribution of the output probabilities $p_{m}({\bf z})$ respects
the condition

$\text{Pr}(p_{m}(\textbf{z})>\delta/N)\geq\gamma,$ (10)

for all possible z. Here $\delta$ and $\gamma$ are positive constants,
$N=2^{L}$ is the dimension of the Hilbert space, and the distribution is
computed over all experimentally implementable unitaries $\\{\hat{U}\\}$.
Equation (10) ensures that most of the basis states have non-zero probability
of being measured. Following the anti-concentration condition, one can show
that if $\hat{U}$ follows the COE, the existence of any efficient classical
sampler $\mathcal{C}$ would imply the collapse of the complexity hierarchy to
the third level [11].

In the case where the dynamics is such that the output state
$|\psi_{m}\rangle$ has an equal probability of being anywhere in the Hilbert
space, its output distribution Pr$(p_{m}({\bf z}))$ follows the PT
distribution [10],

$\text{PT}(p)=Ne^{-Np}\qquad{\rm for}\qquad N\gg 1,$ (11)

and the dynamics satisfies the anti-concentration condition Eq. (10) with
$\delta=1$ and $\gamma=1/e$. As a consequence, the convergence of the output
distribution toward the PT distribution is a key signature of quantum
supremacy [24, 10, 2]. The difference between these two distributions can be
measured by the KLD [39], defined as ($p=p_{m}({\bf z})$ for readability)

$\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))\equiv\sum_{p}\text{Pr}(p)\log\left(\frac{\text{Pr}(p)}{\text{PT}(p)}\right)\geq
0.$ (12)

The KLD is zero only when $\text{Pr}(p)=\text{PT}(p)$ for all $p$.

Figure 3: Reaching the PT distribution. (a)
$\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$ with $p=p_{m}({\bf z})$ as a
function of the driving cycle $m$ for the driven thermalized phase
($W=3J,\omega=8J$; dark red squares) , the driven MBL phase
($W=30J,\omega=8J$; blue circles) and the prethermalization regime
($W=4J,\omega=20J$; green triangles). (b)-(d) Output probability distributions
${\rm Pr}(Np)$ for $m=10$ in the three different regimes. (b) ${\rm Pr}(Np)$
in the driven thermalized phase ($\omega=8J$) is in a good agreement with the
PT distribution (black line) while it is not in (c) for the prethermalization
regime and in (d) for the driven MBL phase. The other parameters are $L=9$
with $N=2^{L}$, $B_{0}=-\delta B=1.25J$ and each data points has been computed
over $100$ disorder realizations.

In Fig. 3, we plot the KLD between the output distribution Pr$(p=p_{m}({\bf
z}))$ after $m$ driving cycles and the PT distribution PT$(p=p_{m}({\bf z}))$
for the three regimes of a driven Ising chain with $L=9$. In panel (a), we
show the evolution of the KLD as a function of $m$ and see that only in the
driven thermalised phase Pr$(p)$ reaches the PT distribution after a number of
cycles that increases as the driving frequency increases.

###  III.2 Computational complexity of the MBL phase and the prethermalization
regime

In the MBL phase, only a small fraction of the Hilbert space near the initial
state is explored. Consequently, given an initial product state,
$|\psi_{0}\rangle=\otimes_{i}^{L}|\uparrow\rangle$, only a limited number of
output bitstrings can be measured after the quantum evolution and the anti-
concentration condition can never be satisfied (see Appendix A for more
details). Consequently, the existence of a classical sampler $\mathcal{C}$
capable of providing bitstrings from a distribution additively closed to
$p_{m}({\bf z})$ is not prohibited by the accepted hierarchy of the complexity
classes.

The more involved question of whether the simulation of the MBL dynamics is
classically intractable or not has already been studied extensively [35, 36,
37, 34] and the answer lies in the properties of the bipartite entanglement
entropy. For an isolated system described by $|\Psi\rangle$, the bipartite
entanglement entropy $S_{e}$ is obtained by dividing the system into two
subsystems $\mathcal{S}$ and $\mathcal{B}$ and calculating the von Neumann
entropy of one of the reduced density matrix, i.e. $S_{e}=-{\rm
Tr}_{\mathcal{S}}\\{\hat{\rho}_{\mathcal{S}}\log_{2}\hat{\rho}_{\mathcal{S}}\\}$
with $\hat{\rho}_{\mathcal{S}}={\rm
Tr}_{\mathcal{B}}\\{|\Psi\rangle\langle\Psi|\\}$. In the MBL phase, all
eigenstates of $\hat{U}$ obey the “area law” [34], meaning that $S_{e}$ is
independent of the subsystem size in 1D systems. This property suggests an
efficient classical representation using matrix product states (MPS) with low
bond dimensions and the use of density matrix renormalization group (DMRG)
techniques [37, 52]. Moreover, given an initial product state, the
entanglement entropy evolves logarithmically with time [53, 54], allowing
efficient and accurate classical simulation of long-time dynamics. Therefore,
further supporting that driven MBL dynamics never converges towards the PT
distribution which would represent infinite-temperature thermalization and
would imply maximal entanglement entropy. In Appendix A, we analytically show
that the MBL dynamics can be efficiently approximated with MPS up to additive
errors.

The scenario in the prethermalized regime where $\hat{U}$ tends towards the
GOE statistics is very different and many questions remain open. For example,
the entanglement entropy of the eigenstates of $\hat{U}_{\rm GOE}$ obeys the
“volume law” and an initial product state would see its entanglement entropy
increase linearly in time [55]. Both properties force an efficient MPS
representation to fail. On the other hand, the fact that the output
distribution does not follow the PT distribution leaves uncertain whether or
not the output distribution anti-concentrates [see Eq. (10)]. The complexity
class associated with sampling from a quantum state in the prethermalized
regime thus remains unclear.

Figure 4: Probing the phase diagram with the
$\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$. (a)
$\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$ with $p=p_{m}({\bf z})$ after
$m=10$ driving cycles as a function of the disorder strength $W/J$ and the
driving frequency $\omega/J$. The solid black line represents the parameters
for which $\langle r\rangle_{\hat{U}}=0.51$ as a comparison with Fig. 2 (a).
(b) A cut along the frequency $\omega=8J$ (vertical dashed line in (a))
comparing the $\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$ (blue dot, left
y-axis) with $\langle r\rangle_{\hat{U}}$ (gray squares, right y-axis) as a
function of $W/J$ along the phase transition from driven thermalized to driven
MBL. (c) A cut along $W=4J$ (horizontal dashed line in (a)) as a function of
$\omega/J$ demonstrates the transition towards the prethermalization regime.
This time, we compare the $\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$ with
$\langle r\rangle_{\hat{U}}$ and $\langle r\rangle_{\hat{U}_{0}}$ (see Fig. 2
(c) and main text). The other parameters are as in Fig. 2.

##  IV Quantum phases from supremacy signatures: A new order parameter

Now that we have presented the different phases and discussed the
computational complexity to sample from their dynamics, we analyse how quantum
supremacy signatures, here taken to be the KLD to PT, can be directly used to
probe the dynamical phase diagram of the driven Ising chain. The motivation is
not only to provide an alternative order parameter, compared to the mean level
spacing $\langle r\rangle$ and the population imbalance for example, but also
to revisit these phase transitions from a computational complexity point of
view.

###  IV.1 Phase diagram

The protocol consists of initializing the system in a product state
$|\psi_{0}\rangle=\otimes_{i}^{L}|\uparrow\rangle$, letting it evolve for $m$
driving cycles and compute its output probability distribution Pr$(p)$
[$p=p_{m}({\bf z})$]. We repeat the procedure for 100 disorder realizations
and average the results. In Fig. 4 (a), we construct a “phase diagram” from
the KLD of the output distributions to PT as a function of the disorder
strength $W/J$ and the drive frequency $\omega/J$. The solid black lines
represent the parameter boundaries with $\langle r\rangle_{\hat{U}}=0.51$ as
in Fig. 2 for easier comparison. From the low-frequency and low-disorder
enclosed region corresponding to $\langle r\rangle_{\hat{U}}>0.51$, we see
that the convergences of Pr$(p)$ towards the PT distribution at $m=10$
coincide with the convergence of the unitary $\hat{U}$ towards the COE. The
divergence from the PT distribution as $W\rightarrow 0$ is again a consequence
of the non interacting nature of the model described by Eq. (1) with $W=0$.
The PT is also not reached in the high-frequency enclosed region where
$\hat{U}$ converges towards the GOE in the prethermalized regime. As discussed
earlier, the computational complexity associated with this prethermalized
regime has yet to be formally identified.

In Fig. 4 (b) and (c), we show a cut along $\omega=8J$ and $W=4J$
respectively. We explicitly see the ability of the KLD to PT to capture the
changes in $\langle r\rangle_{\hat{U}}$ and thus the phase transition from
driven thermalized to MBL (panel b) and to the prethermalized regime (panel
c). Furthermore, it shows that the output probability distribution starts
converging towards PT when the dynamics stop being accurately described by the
time-averaged Hamiltonian $\hat{H}_{0}$.

###  IV.2 Porter-Thomas distribution and infinite temperature thermalization

So far, the discussions concerning the computational complexity of sampling
from the driven thermalized phase mostly came from a complexity theory point
of view (including Ref. [24]) with the only intuition that the PT distribution
is a signature of the dynamics with a capability to explore its entire Hilbert
space. Here we provide additional physical intuitions by making a direct link
between the PT distribution and thermalization to infinite temperature. This
relation highlights the crucial role of the drive in providing the required
energy to reach the PT distribution in our finite-size model.

In an isolated quantum system $\mathcal{S}_{\rm tot}$ described by a pure
state $|\Psi\rangle$, the concept of thermalization is usually understood as
the emergence of homogeneous statistical properties of the reduced density
matrices $\hat{\rho}_{\mathcal{S}}={\rm
Tr}_{\mathcal{B}}|\Psi\rangle\langle\Psi|\equiv\exp{(-\hat{H}_{\mathcal{S}}/k_{B}T_{\rm
eff})}/\mathcal{Z}$. Here, $\mathcal{S}$ denotes any “small” subsystems with
corresponding Hamiltonian $\hat{H}_{\mathcal{S}}$ of dimensions
$N_{\mathcal{S}}$ and partition function $\mathcal{Z}$ such that
$\mathcal{S}_{\rm tot}=\mathcal{S}+\mathcal{B}$ and $N_{\mathcal{B}}\gg
N_{\mathcal{B}}$. The effective temperature $T_{\rm eff}$ is defined from the
reduced density matrices and should become independent of frequency and of the
subsystem choice as the full system reaches thermalization [55, 28]. By
writing
$|\Psi\rangle=\sum_{i}^{N_{\mathcal{S}}}\sum_{j}^{N_{\mathcal{B}}}c_{ij}|i_{\mathcal{S}}\rangle|j_{\mathcal{B}}\rangle$,
where $|j_{\mathcal{S}}\rangle$ and $|i_{\mathcal{B}}\rangle$ are basis states
of the subsystem and the bath respectively with $c_{ij}\in{\mathbb{C}}$, we
obtain

$\hat{\rho}_{\mathcal{S}}=\sum_{i,j=1}^{N_{\mathcal{S}}}\left[\sum_{k=1}^{N_{\mathcal{B}}}c_{ik}c^{*}_{jk}\right]|i_{\mathcal{S}}\rangle\langle
j_{\mathcal{S}}|.$ (13)

In the case where $N_{\mathcal{B}}\gg 1$ and the output probability
distribution of $|\Psi\rangle$ follows the PT distribution given in Eq. (11),
it is straightforward to show that $\langle
i_{\mathcal{S}}|\hat{\rho}_{\mathcal{S}}|j_{\mathcal{S}}\rangle=\sum_{k=1}^{N_{\mathcal{B}}}c_{ik}c^{*}_{jk}=\delta_{ij}/N_{\mathcal{S}}$
with $\delta_{ij}$ being the Kronecker delta function. The vanishing off-
diagonal elements come as a consequence that the amplitudes $c_{ik}$
underlying the PT distribution are random complex numbers subjected to the
normalization constraint $\langle\Psi|\Psi\rangle=1$. This result leads to an
effective infinite temperature $T_{\rm eff}\rightarrow\infty$. Furthermore,
given that $\hat{\rho}_{\mathcal{S}}\sim\hat{1}$ has to hold for all
subsystems with $N_{\mathcal{B}}\gg 1$ to reach thermalization, $c_{ik}$ being
random complex numbers is the only possible structure. In other words, a
quantum state following PT corresponds to an infinite temperature thermalized
state and vice versa.

Interestingly, this link between PT distribution and infinite temperature
thermal phase suggests an alternative quantum supremacy signature that could
be probed to observe the phase transitions. One could consider to use the
difference between the entanglement entropy $S_{e}$ (see section III.2) of all
possible subsystems $\\{\mathcal{S}\\}$ and the corresponding maximal
(infinite temperature) entropy $S_{\rm max}=\log_{2}(N_{\mathcal{S}})$ as an
order parameter. More details concerning the entanglement entropy as a
possible alternative route are given in Appendix B.

###  IV.3 Drive-induced long-range multi-body interactions

More physical intuitions can be gained by highlighting the crucial role of the
effective long-range multi-body interactions induced by the drive in
facilitating a more efficient exploration of the Hilbert space. To see this,
it is insightful to investigate the time-independent Floquet Hamiltonian
$\hat{H}_{F}$. For interacting systems, it is generally impossible to find an
analytic form for $\hat{H}_{F}$ and one needs to resort to approximations. The
most common approach is using the Magnus expansion to expand
$\exp\left[-iH_{F}T\right]$ in a power series of $E_{c}/\omega$, where $E_{c}$
is a characteristic energy of the Hamiltonian $\hat{H}(t)$ [56, 57, 58]. Doing
so, one can write the Floquet Hamiltonian as
$\hat{H}_{F}=\sum_{l=0}^{\infty}\hat{H}_{F}^{(l)}$, where the two first terms,
for example, read

$\displaystyle\hat{H}_{F}^{(0)}=\frac{1}{T}\int_{0}^{T}d\tau_{1}\hat{H}(\tau_{1})=\hat{H}_{0},$
(14)
$\displaystyle\hat{H}_{F}^{(1)}=\frac{1}{2iT}\int_{0}^{T}d\tau_{1}\int_{0}^{\tau_{1}}d\tau_{2}[\hat{H}(\tau_{1}),\hat{H}(\tau_{2})].$
(15)

The first term is the time average of $\hat{H}(t)$ which equals $\hat{H}_{0}$
and the next correction $H_{F}^{(2)}$ is computed explicitly in Appendix C. In
the limit of infinite frequencies $\omega\rightarrow\infty$, $\hat{H}_{F}$
converges to $\hat{H}_{0}$, while the series diverges in the low frequency
regime $\omega<E_{c}$.

The crucial point is that higher-order contributions generally include multi-
body and longer-range effective interactions. For example, $\hat{H}_{F}^{(2)}$
includes three-body interaction terms of the form
$\hat{H}_{i}\sim\frac{J^{2}}{\omega^{2}}\delta
B\hat{\sigma}^{z}_{i-1}\hat{\sigma}^{x}_{i}\hat{\sigma}^{z}_{i+1}$ [cf. Eq.
(23) in Appendix C]. This tendency suggests that in the limit where the series
diverges, the dynamics is governed by infinitely-long range multi-body
interactions. In Fig. 2 (b) we show that as the driving frequency decreases,
the level statistics of $\hat{U}$ starts converging toward the COE as the
higher-order terms in the Magnus expansion start playing an increasingly
important role. This is captured by comparing the statistics of the full
evolution $\hat{U}$ with the statistics of the evolution $\hat{U}_{0}$ solely
generated by the averaged Hamiltonian $\hat{H}_{0}$. The same is observed in
Fig. 4 (a) where the output probability distribution Pr$(p)$ starts converging
toward the PT distribution only once the evolution under $\hat{H}_{0}$ fails
at approximating the exact evolution. This observation sheds light on the
advantages of periodically driven systems to explore larger regions of their
Hilbert space by highlighting the effectively enhanced connectivity of the
model, which leads to faster growth of entanglement. Note that two-body long-
range interactions in 1D systems have been shown to slow down the entanglement
entropy growth compared to nearest-neighbour interacting models [59]. This
results further highlights the importance of the multi-body nature of the
induced interactions. In Appendix D, we discuss how consequently driven analog
quantum simulators should reach the quantum supremacy regime more rapidly
compared to digital random quantum circuits with the same topology.

##  V Conclusion and outlooks

In this work, we have described how key signatures of quantum supremacy in
driven analog interacting systems can be directly used to probe its dynamical
phase diagram. More precisely, we have focused on showing how the KLD between
the distribution of the output state probabilities and the PT distribution can
play the role of an order parameter capturing the transitions from the driven
thermalized phase to the driven MBL phase and the Floquet prethermalized
regime. In this context, we have highlighted the clear connection between the
effective multi-body long-range interactions generated by the drive and the
capability of the system to reach the PT distribution. Compared to the
commonly measured population imbalance [34, 42, 43, 44] which captures the
diffusion of an initial state in real space, the KLD to PT tells us about the
diffusion over the Hilbert space. Beside proposing an alternative order
parameter for dynamical phase transitions, this work adopts a different angle
more oriented towards the complexity theory. It represents a direct
application of the very same protocol proposed to demonstrate quantum
supremacy in generic driven many-body systems. An interesting next step would
be to analyse in more details the computational complexity aspects of the
Floquet prethermalized regime.

##  VI Acknowledgement

This research is supported by the National Research Foundation, Prime
Minister’s Office, Singapore and the Ministry of Education, Singapore under
the Research Centres of Excellence programme. It was also partially funded by
Polisimulator project co-financed by Greece and the EU Regional Development
Fund, Ninnat Dangniam is supported by the National Natural Science Foundation
of China (Grant No. 11875110).

##  Appendix A MBL phase

We here provide a more detailed mathematical description of the MBL phase. We
focus on demonstrating that the probability distribution $p_{m}({\bf z})$ of a
quantum state evolving under its dynamics failed to satisfy the anti-
concentration condition given in Eq. (10). Finally, we discuss the feasibility
to efficiently construct a classical sampler $\mathcal{C}$ capable of
providing bitstrings from a distribution $q({\bf z})$ that is additively
closed to $p_{m}({\bf z})$ as defined in Eq. (9).

###  A.1 General framework

The standard approach to describe the MBL phase is via a complete set of
quasi-local integral of motions [60, 61, 62, 34], defined as
$[\hat{\tau}^{z}_{i},\hat{H}_{\rm MBL}]=0$, where

$\hat{\tau}^{z}_{i}=Z\hat{\sigma}^{z}_{i}+\sum_{m=1}e^{-m/\xi}\hat{O}^{(m)}_{i}.$
(16)

Here, $Z\approx 1$, $\xi\ll 1$ is the characteristic localization length and
$\hat{O}^{(m)}_{i}$ are generic operators describing interaction between up to
$(2m+1)-$body over sites from $i-m$ to $i+m$. A generic MBL Hamiltonian then
adopts the form

$\displaystyle\hat{H}_{\rm MBL}=$
$\displaystyle\sum_{i}h_{i}\hat{\tau}^{z}_{i}+\sum_{i<j}J_{0}e^{-|i-j|/\xi^{\prime}}\hat{\tau}^{z}_{i}\hat{\tau}^{z}_{j}$
$\displaystyle+\sum_{i<j<k}J_{0}e^{-|i-k|/\xi^{\prime}}\hat{\tau}^{z}_{i}\hat{\tau}^{z}_{j}\hat{\tau}^{z}_{k}+...,$
(17)

where $h_{i}\in[0,W]$ are on-site energies uniformly drawn from a random
distribution, $J_{0}$ is an effective interaction strength and
$\xi^{\prime}\sim\xi$ is an additional localization length. For periodically
driven systems, the standard approach is to cast the Floquet Hamiltonian
$\hat{H}_{F}$ onto the form of Eq. (A.1) [33].

###  A.2 Concentration of output distribution

By inserting Eq. (16) into Eq. (A.1), The MBL Hamiltonian can be expanded as a
perturbation series

$\displaystyle\hat{H}_{\rm
MBL}=\hat{H}_{0}+\sum_{m=1}^{\infty}\sum_{i=1}^{L}\epsilon^{m}\hat{V}^{(m)}_{i},$
(18)

where $\epsilon=e^{-1/\xi}\ll 1$ (we took $\xi^{\prime}=\xi$ for simplicity),
$\hat{H}_{0}=\sum_{i}h_{i}Z\hat{\sigma}^{z}_{i}$ is the non-interacting part
and $\hat{V}^{(m)}_{i}$ are local perturbations acting on at most sites from
$i-m$ to $i+m$. At first order perturbation ($m=1$), $\hat{V}^{(1)}_{i}$ can
flip at max 3 neighbouring spins and can thus only mix an unperturbed product
state of $L$ spins $|{\bf z}_{0}\rangle$ with up to $4L$ different basis
states (neglecting the boundary corrections). The same logic applies at second
order perturbation where this time $|{\bf z}_{0}\rangle$ becomes a linear
superposition of up to $16L$ additional basis states. As the order of
perturbation increases, a polynomial number (poly($L$)) of additional basis
states have to be taken into account. This is a consequence of the local
nature of each perturbation term. In contrast, extended interactions, as it
would be the case in the thermal phase, would lead to an exponential increase
of additional basis states. Assuming that an initial product state evolves
deep inside the MBL phase where only few orders of perturbation is required,
the output state would only have a polynomial number of basis state populated.
Since the anti-concentration condition in Eq. (10) requires most of the
$2^{L}$ basis states to be populated, the output distribution in the MBL phase
can not satisfy it.

Figure 5: Probing phase diagram with entanglement entropy. (a) The average
entanglement entropy over six different subsystems, each containing three
sites chosen randomly, is used as an alternative quantum supremacy signature
to revisit the phase diagram. The vertical and horizontal cuts show the
average entanglement entropy as a function of (b) $W/J$ and (c) $\omega/J$.
The shaded orange region indicates one standard deviation. The horizontal
dashed lines represent the maximal entropy at the infinite temperature $S_{\rm
max}=S_{e}(T_{\rm eff}=\infty)$. The KLD to PT
$\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$ (blue circles) and the average
level spacing $\langle r\rangle_{\hat{U}}$ (grey squares) are also plotted for
comparison. All parameters are identical to those in Fig. 4. The low variance
in the thermalized phase indicates the homogeneity.

###  A.3 Complexity of sampling task

Previous works have demonstrated the ability of MPS states to efficiently
approximate all eigenstates of any MBL Hamiltonians [35, 34]. In what follows
we denote an exact eigenstate of $\hat{H}_{\rm MBL}$ with eigenenergy $E_{k}$
as $|\varphi_{k}\rangle$ and the associated MPS approximation as $|{\rm
MPS}_{k}\rangle$. We define the error made by the approximation, $\nu\ll 1$,
as

$\displaystyle|{\rm MPS}_{k}\rangle$
$\displaystyle=(1-\nu)|\varphi_{k}\rangle+\sqrt{2\nu-\nu^{2}}\sum_{k^{\prime}\neq
k}a_{kk^{\prime}}|\varphi_{k^{\prime}}\rangle,$
$\displaystyle\approx|\varphi_{k}\rangle+\sqrt{2\nu}\sum_{k^{\prime}\neq
k}a_{kk^{\prime}}|\varphi_{k^{\prime}}\rangle,$ (19)

where $\sum_{k^{\prime}\neq k}a_{kk^{\prime}}$ captures the overlap over all
other eigenstates $|\varphi_{k^{\prime}}\rangle$ with $\sum_{k^{\prime}\neq
k}|a_{kk^{\prime}}|^{2}=1$. The mean energy associated with $|{\rm
MPS}_{k}\rangle$ can be written as

$\displaystyle\tilde{E}_{k}$ $\displaystyle=\langle{\rm MPS}_{k}|\hat{H}_{\rm
MBL}|{\rm MPS}_{k}\rangle,$ $\displaystyle\approx
E_{k}+2\nu\sum_{k^{\prime}\neq k}|a_{kk^{\prime}}|^{2}E_{k^{\prime}},$
$\displaystyle\equiv E_{k}+2\nu\Lambda_{k},$ (20)

where $\Lambda_{k^{\prime}}$ represents the average over the entire energy
spectrum following the “error distribution” $\sum_{k^{\prime}\neq
k}|a_{kk^{\prime}}|^{2}$. The bond dimension of the MPS states required to
maintain a fix level of precision $\nu$ have been shown to scale polynomially
with the system size $L$ [35, 34]; the exact scaling depends on the specifics
of $\hat{H}_{\rm MBL}$. The approximated output probability distribution
$q_{t}(\textbf{z})$, after evolving from the initial state $|{\bf
z}_{0}\rangle$ for a time $t$, then reads

$\displaystyle q_{t}(\textbf{z})$
$\displaystyle=\left|\langle\textbf{z}|\left(\sum_{k}e^{-i\tilde{E}_{k}t}|\rm{MPS}_{k}\rangle\langle\rm{MPS}_{k}|\right)|\textbf{z}_{0}\rangle\right|^{2},$
(21)
$\displaystyle\approx\left|\langle\textbf{z}|\sum_{k}e^{-i(E_{k}+2\nu\Lambda_{k})t}\left(|\varphi_{k}\rangle\langle\varphi_{k}|+\sqrt{2\nu}\sum_{k^{\prime}\neq
k}\left[a_{kk^{\prime}}|\varphi_{k^{\prime}}\rangle\langle\varphi_{k}|+{\rm
H.c.}\right]\right)|\textbf{z}_{0}\rangle\right|^{2},$
$\displaystyle\approx\left|\langle\textbf{z}|\sum_{k}e^{-iE_{k}t}\left([1-2i\nu\Lambda_{k}t]|\varphi_{k}\rangle\langle\varphi_{k}|+\sqrt{2\nu}\sum_{k^{\prime}\neq
k}\left[a_{kk^{\prime}}|\varphi_{k^{\prime}}\rangle\langle\varphi_{k}|+{\rm
H.c.}\right]\right)|\textbf{z}_{0}\rangle\right|^{2},$ $\displaystyle\approx
p(\textbf{z})+\left[\sqrt{2\nu}\sum_{k,k^{\prime}}e^{-i(E_{k}-E_{k^{\prime}})t}\langle\textbf{z}|\hat{P}_{k}|\textbf{z}_{0}\rangle\langle\textbf{z}_{0}|\hat{M}_{k^{\prime}}|\textbf{z}\rangle+2i\nu
t\sum_{k,k^{\prime}}e^{-i(E_{k}-E_{k^{\prime}})t}\Lambda_{k^{\prime}}\langle\textbf{z}|\hat{P}_{k}|\textbf{z}_{0}\rangle\langle\textbf{z}_{0}|\hat{P}_{k^{\prime}}|\textbf{z}\rangle+c.c.\right].$

Here $\hat{P}_{k}=|\varphi_{k}\rangle\langle\varphi_{k}|$,
$\hat{M}_{k}=\sum_{k^{\prime}\neq
k}a_{kk^{\prime}}|\varphi_{k^{\prime}}\rangle\langle\varphi_{k}|+{\rm H.c.}$
and the exact distribution is given by
$p_{t}(\textbf{z})=\left|\sum_{k}e^{-iE_{k}t}\langle\textbf{z}|\hat{P}_{k}|\textbf{z}_{0}\rangle\right|^{2}$.
In the above expression, we keep only the leading order in $\nu$ and $\nu t$
and we considered the limit $\nu\Lambda_{k}t\ll 1$.

In the short time limit where the last term in the last line of Eq. (21) can
be neglected, the error made in approximating $p_{t}(\textbf{z})$ is
proportional to
$\sqrt{2\nu}\sum_{k}e^{-iE_{k}t}\langle\textbf{z}|\hat{P}_{k}|\textbf{z}_{0}\rangle$.
Following the discussion in the previous subsection, we know that for a
dynamics deep in the MBL phase, only a small fraction [poly$(L)/2^{L}$] of
$\langle\textbf{z}|\hat{P}_{k}|\textbf{z}_{0}\rangle$ will be nonzero and at
max proportional to $e^{-1/\xi}$. Therefore, the additive error condition in
Eq. (9), i.e.
$\sum_{\textbf{z}}|p_{t}(\textbf{z})-q_{t}(\textbf{z})|\leq\beta_{0}$, is
expected to be reachable with arbitrary small $\beta_{0}$ for bond dimensions
of the MPS states that scale polynomially with the system size. For longer
evolution times, the term proportional to $\nu t$ will start to become
increasingly important. It will inevitably reach a point where any small
initial mistake $2\nu\Lambda_{k}$ made estimating the eigenenergies
accumulates, leading to a large error in the phases. While any approximation
method should fail in this very large time limit, we note that the same is
bounded to happen for a quantum simulator; any small amount of dephasing or
decay will inevitably accumulate beyond the required accuracy in the large
time limit.

##  Appendix B Entanglement entropy

As discussed in section IV.2, dividing a closed system into two parts,
${\mathcal{S}}$ and ${\mathcal{B}}$ of size $N_{\mathcal{S}}$ and
$N_{\mathcal{B}}$ respectively, and having a reduced density matrix
$\hat{\rho}_{\mathcal{S}}=\hat{1}/N_{\mathcal{S}}$ for all possible subsystems
$\mathcal{S}$ corresponds to an infinite-temperature thermalized state. Here
the reduced density matrix is $\hat{\rho}_{\mathcal{S}}={\rm
Tr}_{\mathcal{B}}|\Psi\rangle\langle\Psi|$ where $|\Psi\rangle$ is the state
of the full closed system. In that case, the output probability distribution
of the quantum state $|\Psi\rangle$ follows the PT distribution and
consequently fulfil the anti-concentration condition given in Eq. (10). This
means that an alternative route to demonstrate the anti-concentration
condition could be to measure the reduced density matrix of all subsystems
$\mathcal{S}$ and show that $\hat{\rho}_{\mathcal{S}}=\hat{1}/N_{\mathcal{S}}$
homogeneously over the whole system. We note that entanglement entropy of
undriven systems has been experimentally measured in different quantum
platforms [34, 55, 44, 63].

In Fig. 5, we plot the entanglement entropy $S_{e}={\rm
Tr}\\{\hat{\rho}_{\mathcal{S}}\log\hat{\rho}_{\mathcal{S}}\\}$ as the function
of $\omega$ and $W$ averaged over multiple subsystems. We choose six different
subsystems consisting of three sites randomly selected and apply the same
procedure as in the main text. We see that the system approaches the maximal
entropy associated with the infinite temperature state, $S_{\rm
max}=\log_{2}(N_{\mathcal{S}})$ when the output probability follows PT and the
statistics of $\hat{U}$ converges toward the COE. The homogeneity of the
thermalized phase is indicated by the low variance of $S_{e}$ over the
different subsystems.

##  Appendix C Magnus expansion of $\hat{H}(t)$.

We here explicitly compute the Magnus expansion up to the second order
correction. The first term reads

$\displaystyle\hat{H}_{F}^{(1)}$
$\displaystyle=\frac{1}{2iT}\int_{t_{0}}^{T+t_{0}}d\tau_{1}\int_{t_{0}}^{\tau_{1}}d\tau_{2}\left[\hat{H}(\tau_{1}),\hat{H}(\tau_{2})\right]$
(22)
$\displaystyle=\frac{1}{2iT}\int_{t_{0}}^{T+t_{0}}d\tau_{1}\int_{t_{0}}^{\tau_{1}}d\tau_{2}\left[\hat{H}_{0},\hat{H}_{d}(\tau_{2})-\hat{H}_{d}(\tau_{1})\right]$
$\displaystyle=\frac{2i\delta B\sin{(\omega
t_{0}})}{\omega}\left[\sum_{j=1}^{L}h_{j}\hat{\sigma}^{y}_{j}+J\sum_{j=1}^{L-1}\left(\hat{\sigma}^{y}_{j}\hat{\sigma}^{z}_{j+1}+\hat{\sigma}^{z}_{j}\hat{\sigma}^{y}_{j+1}\right)\right].$

The choice of $t_{0}=0$ is arbitrary and can lead to different forms of the
Floquet Hamiltonian $\hat{H}_{F}$. Consistently with the main text, we fix
$t_{0}=0$, leading to $\hat{H}_{F}^{(1)}=0$. The second order correction term
adopts the form

$\displaystyle\hat{H}_{F}^{(2)}=$
$\displaystyle-\frac{1}{6T}\int_{t_{0}}^{T+t_{0}}d\tau_{1}\int_{t_{0}}^{\tau_{1}}d\tau_{2}\int_{t_{0}}^{\tau_{2}}d\tau_{3}\times\left(\left[\hat{H}(\tau_{1}),\left[\hat{H}(\tau_{2}),\hat{H}(\tau_{3})\right]\right]+(1\Leftrightarrow
3)\right),$ $\displaystyle=$ $\displaystyle\frac{-4\delta
B}{\omega^{2}}\left\\{\sum_{j=1}^{L}h_{j}^{2}\hat{\sigma}^{x}_{j}+2J\sum_{j=1}^{L-1}\left(h_{j}\hat{\sigma}^{x}_{j}\hat{\sigma}^{z}_{j+1}+h_{j+1}\hat{\sigma}^{z}_{j}\hat{\sigma}^{x}_{j+1}\right)+J^{2}\sum_{j=1}^{L-1}\left(\hat{\sigma}^{x}_{j}+\hat{\sigma}^{x}_{j+1}\right)+2J^{2}\sum_{j=1}^{L-2}\left(\hat{\sigma}^{z}_{j}\hat{\sigma}^{x}_{j+1}\hat{\sigma}^{z}_{j+2}\right)\right\\}$
$\displaystyle+\frac{4B_{0}\delta
B}{\omega^{2}}\left\\{\sum_{j=1}^{L}h_{j}\hat{\sigma}^{z}_{j}+2J\sum_{j=1}^{L-1}\left(\hat{\sigma}^{z}_{j}\hat{\sigma}^{z}_{j+1}-\hat{\sigma}^{y}_{j}\hat{\sigma}^{y}_{j+1}\right)\right\\}.$
(23)

As discussed in the main text, higher-order correction includes longer-range
interactions involving increasing number of spins with coupling strength that
scale proportionally to $E_{c}/\omega$. Multiple works have been dedicated to
the study of the Magnus expansion convergence in driven systems by estimating
the relevant energy scale $E_{c}$ [64, 65].

Figure 6: Benchmarking quantum supremacy with analog and digital approaches:
(a) a circuit diagram of the analog driven quantum Ising chain with $m$
driving cycles. (b) a circuit diagram of random quantum circuits containing
$m$ layers of random single qubits in the set $\\{\sqrt{X},\sqrt{Y},T\\}$ and
controlled-Z gates. (c) ((d)) $\text{KL}(\text{Pr}(p)\parallel\text{PT}(p))$
as a function of $m$ and the size of the system $L$ for the analog (digital)
approach. (e) The output distribution, weighted by $N=2^{L}$, in the analog
case with $L=11$ and $m=40$. The yellow line is the exact PT distribution. (f)
KL divergence at the long time limit of the analog and the digital cases as a
function of $L$. ($W=3J$, $B_{0}=-\delta B=1.25J$, $\omega=8J$ and 500
disorder realizations.)

##  Appendix D Comparison between analog and digital approaches in reaching
quantum supremacy

In this appendix, we compare the time it takes for a quantum state to evolve
toward the PT distribution under the dynamics of the driven analog Ising spin
chain and the 1D digital random quantum circuits [10]. To keep the comparison
fair, we consider the same interaction topology, i.e. a 1D system with
nearest-neighbour couplings and open boundary conditions. We show that the
effective long-range multi-body interactions induced by the drive allow the
analog quantum system to converge considerably faster in the regime where
quantum supremacy can be demonstrated.

The model we use to simulate the random quantum circuits is based on Ref.[10]
and consists of $m$ layers of gates acting on $L$ qubits, as depicted in Fig.
6(b). Each layer consists of a series of single-qubit gates followed by a
series of controlled-Z gates. Moreover, a layer of Hadamard gates is added
once at the very beginning of the circuit. Each single-qubit gates are chosen
randomly from the set $\\{\sqrt{X},\sqrt{Y},T\\}$, where $\sqrt{X}$
($\sqrt{Y}$) represents a $\pi/2$ rotation around the $X$ ($Y$) axis of the
Bloch sphere and $T$ is a non-Clifford gate represented by the diagonal matrix
$\\{1,e^{i\pi/4}\\}$. As in the analog case, we measure the output
distribution ${\rm Pr}(p_{m}({\bf z}))$ averaged over $D$ distinct
realizations of this model circuit. In order to compare the evolution time, we
set the driving period of the analog Ising chain to be approximately equal to
the duration of a single layer of gates ($m=1$) of the digital circuit. Since
single-qubit gates can usually be implemented much faster compared to two-
qubit gates, we assume that each layer takes the same time as a controlled-Z
gate from a standard $J\sigma^{z}_{i}\sigma^{z}_{i+1}$ coupling, i.e.
$t_{m}\sim\pi/4J$. This results in a driving frequency $\omega=8J$.

In Fig. 6(c)-(d), we plot the $\text{KLD}(\text{Pr}(p)\parallel\text{PT}(p))$
with $p=p_{m}({\bf z})$ for the driven Ising chain and the digital random
circuits, respectively, as a function of the driving cycle $m$ for different
system sizes $L$. The long time limit is shown in Fig. 6 (e). We observe that,
despite having the same topology, the driven analog system converges to the PT
distribution considerably faster as a results of the enhanced connectivity
induced by the drive. Given the limited coherence times in the current state-
of-the-art quantum platforms, such gain makes the analog approaches more
favorable as testbeds for quantum supremacy. Moreover, the analog approach
does not require precise local control to implement the gates.

## References

  * [1] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando G. S. L. Brandao, David A. Buell, Brian Burkett, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins, William Courtney, Andrew Dunsworth, Edward Farhi, Brooks Foxen, Austin Fowler, Craig Gidney, Marissa Giustina, Rob Graff, Keith Guerin, Steve Habegger, Matthew P. Harrigan, Michael J. Hartmann, Alan Ho, Markus Hoffmann, Trent Huang, Travis S. Humble, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Paul V. Klimov, Sergey Knysh, Alexander Korotkov, Fedor Kostritsa, David Landhuis, Mike Lindmark, Erik Lucero, Dmitry Lyakh, Salvatore Mandrà, Jarrod R. McClean, Matthew McEwen, Anthony Megrant, Xiao Mi, Kristel Michielsen, Masoud Mohseni, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Murphy Yuezhen Niu, Eric Ostby, Andre Petukhov, John C. Platt, Chris Quintana, Eleanor G. Rieffel, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Kevin J. Sung, Matthew D. Trevithick, Amit Vainsencher, Benjamin Villalonga, Theodore White, Z. Jamie Yao, Ping Yeh, Adam Zalcman, Hartmut Neven, and John M. Martinis.  Quantum supremacy using a programmable superconducting processor.  Nature, 574(7779):505–510, 2019.
  * [2] C. Neill, P. Roushan, K. Kechedzhi, S. Boixo, S. V. Isakov, V. Smelyanskiy, A. Megrant, B. Chiaro, A. Dunsworth, K. Arya, R. Barends, B. Burkett, Y. Chen, Z. Chen, A. Fowler, B. Foxen, M. Giustina, R. Graff, E. Jeffrey, T. Huang, J. Kelly, P. Klimov, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, H. Neven, and J. M. Martinis.  A blueprint for demonstrating quantum supremacy with superconducting qubits.  Science, 360(6385):195–199, 2018.
  * [3] Hui Wang, Jian Qin, Xing Ding, Ming-Cheng Chen, Si Chen, Xiang You, Yu-Ming He, Xiao Jiang, L. You, Z. Wang, C. Schneider, Jelmer J. Renema, Sven Höfling, Chao-Yang Lu, and Jian-Wei Pan.  Boson sampling with 20 input photons and a 60-mode interferometer in a $1{0}^{14}$-dimensional hilbert space.  Phys. Rev. Lett., 123:250503, Dec 2019.
  * [4] Han-Sen Zhong, Hui Wang, Yu-Hao Deng, Ming-Cheng Chen, Li-Chao Peng, Yi-Han Luo, Jian Qin, Dian Wu, Xing Ding, Yi Hu, Peng Hu, Xiao-Yan Yang, Wei-Jun Zhang, Hao Li, Yuxuan Li, Xiao Jiang, Lin Gan, Guangwen Yang, Lixing You, Zhen Wang, Li Li, Nai-Le Liu, Chao-Yang Lu, and Jian-Wei Pan.  Quantum computational advantage using photons.  Science, 2020.
  * [5] Daniel J. Brod, Ernesto F. Galvão, Andrea Crespi, Roberto Osellame, Nicolò Spagnolo, and Fabio Sciarrino.  Photonic implementation of boson sampling: a review.  Advanced Photonics, 1(3):1 – 14, 2019.
  * [6] John Preskill.  Quantum Computing in the NISQ era and beyond.  Quantum, 2:79, August 2018.
  * [7] A. W. Harrow and A. Montanaro.  Quantum computational supremacy.  Nature, 549:203, 09 2017.
  * [8] A. P. Lund, Michael J. Bremner, and T. C. Ralph.  Quantum sampling problems, bosonsampling and quantum supremacy.  npj Quantum Information, 3(1):15, 2017.
  * [9] Michael J. Bremner, Richard Jozsa, and Dan J. Shepherd.  Classical simulation of commuting quantum computations implies collapse of the polynomial hierarchy.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467(2126):459–472, 2011.
  * [10] 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.  Nature Physics, 14(6):595–600, 2018.
  * [11] Adam Bouland, Bill Fefferman, Chinmay Nirkhe, and Umesh Vazirani.  On the complexity and verification of quantum random circuit sampling.  Nature Physics, 15(2):159–163, 2019.
  * [12] Ramis Movassagh.  Quantum supremacy and random circuits.  arXiv e-prints, page arXiv:1909.06210, September 2019.
  * [13] S. Aaronson and A. Arkhipov.  The computational complexity of linear optics.  Proceedings of the 43rd annual ACM Symposium on Theory of Computing, (STOC ’11), pages 333–342, 2011.
  * [14] Aaronson Scott.  A linear-optical proof that the permanent is np-hard.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467:3393, 2011.
  * [15] Juan Bermejo-Vega, Dominik Hangleiter, Martin Schwarz, Robert Raussendorf, and Jens Eisert.  Architectures for quantum simulation showing a quantum speedup.  Phys. Rev. X, 8:021010, Apr 2018.
  * [16] Xun Gao, Sheng-Tao Wang, and L.-M. Duan.  Quantum supremacy for simulating a translation-invariant ising spin model.  Phys. Rev. Lett., 118:040502, Jan 2017.
  * [17] Brian Coyle, Daniel Mills, Vincent Danos, and Elham Kashefi.  The born supremacy: quantum advantage and training of an ising born machine.  npj Quantum Information, 6(1):60, 2020.
  * [18] Ramis Movassagh.  Quantum supremacy and random circuits, 2020.
  * [19] Jarrod R. McClean, Sergio Boixo, Vadim N. Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  Nature Communications, 9(1):4812, 2018.
  * [20] Yudong Cao, Jonathan Romero, Jonathan P. Olson, Matthias Degroote, Peter D. Johnson, Mária Kieferová, Ian D. Kivlichan, Tim Menke, Borja Peropadre, Nicolas P. D. Sawaya, Sukin Sim, Libor Veis, and Alán Aspuru-Guzik.  Quantum chemistry in the age of quantum computing.  Chemical Reviews, 119(19):10856–10915, 10 2019.
  * [21] Sam McArdle, Suguru Endo, Alán Aspuru-Guzik, Simon C. Benjamin, and Xiao Yuan.  Quantum computational chemistry.  Rev. Mod. Phys., 92:015003, Mar 2020.
  * [22] Bela Bauer, Sergey Bravyi, Mario Motta, and Garnet Kin-Lic Chan.  Quantum algorithms for quantum chemistry and quantum materials science, 2020.
  * [23] Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.  Parameterized quantum circuits as machine learning models.  Quantum Science and Technology, 4(4):043001, nov 2019.
  * [24] Jirawat Tangpanitanon, Supanut Thanasilp, Marc-Antoine Lemonde, Ninnat Dangniam, and Dimitris G. Angelakis.  Quantum supremacy in driven quantum many-body systems.  arXiv e-prints, page arXiv:2002.11946, February 2020.
  * [25] André Eckardt.  Colloquium: Atomic quantum gases in periodically driven optical lattices.  Rev. Mod. Phys., 89:011004, Mar 2017.
  * [26] C. Monroe, W. C. Campbell, L. M. Duan, Z. X. Gong, A. V. Gorshkov, P. Hess, R. Islam, K. Kim, N. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao.  Programmable Quantum Simulations of Spin Systems with Trapped Ions.  arXiv e-prints, page arXiv:1912.07845, December 2019.
  * [27] L. D’Alessio and M. Rigol.  Long-time behavior of isolated periodically driven interacting lattice systems.  Phys. Rev. X, 4:041048, Dec 2014.
  * [28] Takashi Mori, Tatsuhiko N Ikeda, Eriko Kaminishi, and Masahito Ueda.  Thermalization and prethermalization in isolated quantum systems: a theoretical overview.  Journal of Physics B: Atomic, Molecular and Optical Physics, 51(11):112001, may 2018.
  * [29] Pedro Ponte, Anushya Chandran, Z. Papić, and Dmitry A. Abanin.  Periodically driven ergodic and many-body localized quantum systems.  Annals of Physics, 353:196 – 204, 2015.
  * [30] L. D’Alessio and A. Polkovnikov.  Many-body energy localization transition in periodically driven systems.  Annals of Physics, 333:19 – 33, 2013.
  * [31] P. Ponte, Z. Papić, F. Huveneers, and D. A. Abanin.  Many-body localization in periodically driven systems.  Phys. Rev. Lett., 114:140401, Apr 2015.
  * [32] Achilleas Lazarides, Arnab Das, and Roderich Moessner.  Fate of many-body localization under periodic driving.  Phys. Rev. Lett., 115:030402, Jul 2015.
  * [33] D. A. Abanin, W. D. Roeck, and F. Huveneers.  Theory of many-body localization in periodically driven systems.  Annals of Physics, 372:1 – 11, 2016.
  * [34] Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, and Maksym Serbyn.  Colloquium: Many-body localization, thermalization, and entanglement.  Rev. Mod. Phys., 91:021001, May 2019.
  * [35] M. Friesdorf, A. H. Werner, W. Brown, V. B. Scholz, and J. Eisert.  Many-body localization implies that eigenvectors are matrix-product states.  Phys. Rev. Lett., 114:170505, May 2015.
  * [36] A. Chandran, J. Carrasquilla, I. H. Kim, D. A. Abanin, and G. Vidal.  Spectral tensor networks for many-body localization.  Phys. Rev. B, 92:024201, Jul 2015.
  * [37] Thorsten B. Wahl, Arijeet Pal, and Steven H. Simon.  Efficient representation of fully many-body localized systems using tensor networks.  Phys. Rev. X, 7:021018, May 2017.
  * [38] Jirawat Tangpanitanon, Supanut Thanasilp, Ninnat Dangniam, Marc-Antoine Lemonde, and Dimitris G. Angelakis.  Expressibility and trainability of parameterized analog quantum systems for machine learning applications.  arXiv: Quantum Physics, 2020.
  * [39] S. Kullback and R. A. Leibler.  On information and sufficiency.  Ann. Math. Statist., 22(1):79–86, 03 1951.
  * [40] C. E. Porter and R. G. Thomas.  Fluctuations of nuclear reaction widths.  Phys. Rev., 104:483–491, Oct 1956.
  * [41] This is similar to the growth of the bipartite entanglement entropy in closed systems which also captures the increase in the capacity of the system to explore its Hilbert space [45]. While a large entanglement entropy implies a classically intractable matrix product state (MPS) representation [66], a vanishing KLD from the PT distribution has been shown to be a key condition to prove that any classical approach is intractable [10, 24].
  * [42] Antonio Rubio-Abadal, Jae-yoon Choi, Johannes Zeiher, Simon Hollerith, Jun Rui, Immanuel Bloch, and Christian Gross.  Many-body delocalization in the presence of a quantum bath.  Phys. Rev. X, 9:041014, Oct 2019.
  * [43] Thomas Kohlert, Sebastian Scherg, Xiao Li, Henrik P. Lüschen, Sankar Das Sarma, Immanuel Bloch, and Monika Aidelsburger.  Observation of many-body localization in a one-dimensional system with a single-particle mobility edge.  Phys. Rev. Lett., 122:170403, May 2019.
  * [44] B. Chiaro, C. Neill, A. Bohrdt, M. Filippone, F. Arute, K. Arya, R. Babbush, D. Bacon, J. Bardin, R. Barends, S. Boixo, D. Buell, B. Burkett, Y. Chen, Z. Chen, R. Collins, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen, C. Gidney, M. Giustina, M. Harrigan, T. Huang, S. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, E. Lucero, J. McClean, X. Mi, A. Megrant, M. Mohseni, J. Mutus, M. McEwen, O. Naaman, M. Neeley, M. Niu, A. Petukhov, C. Quintana, N. Rubin, D. Sank, K. Satzinger, A. Vainsencher, T. White, Z. Yao, P. Yeh, A. Zalcman, V. Smelyanskiy, H. Neven, S. Gopalakrishnan, D. Abanin, M. Knap, J. Martinis, and P. Roushan.  Direct measurement of non-local interactions in the many-body localized phase, 2019.
  * [45] Luca D’Alessio, Yariv Kafri, Anatoli Polkovnikov, and Marcos Rigol.  From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics.  Advances in Physics, 65(3):239–362, 2016.
  * [46] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis.  Spectroscopic signatures of localization with interacting photons in superconducting qubits.  Science, 358(6367):1175–1179, 2017.
  * [47] Antonio Rubio-Abadal, Matteo Ippoliti, Simon Hollerith, David Wei, Jun Rui, S. L. Sondhi, Vedika Khemani, Christian Gross, and Immanuel Bloch.  Floquet prethermalization in a bose-hubbard system.  Phys. Rev. X, 10:021044, May 2020.
  * [48] Pranjal Bordia, Henrik Lüschen, Ulrich Schneider, Michael Knap, and Immanuel Bloch.  Periodically driving a many-body localized quantum system.  Nature Physics, 13:460 EP –, 01 2017.
  * [49] Fritz Haake.  Quantum Signatures of Chaos.  Springer International Publishing, US, 2010.
  * [50] Angelo Russomanno, Giuseppe E Santoro, and Rosario Fazio.  Entanglement entropy in a periodically driven ising chain.  Journal of Statistical Mechanics: Theory and Experiment, 2016(7):073101, jul 2016.
  * [51] Dominik Hangleiter, Juan Bermejo-Vega, Martin Schwarz, and Jens Eisert.  Anticoncentration theorems for schemes showing a quantum speedup.  Quantum, 2:65, May 2018.
  * [52] Ulrich Schollwöck.  The density-matrix renormalization group in the age of matrix product states.  Annals of Physics, 326(1):96 – 192, 2011.  January 2011 Special Issue.
  * [53] Jens H. Bardarson, Frank Pollmann, and Joel E. Moore.  Unbounded growth of entanglement in models of many-body localization.  Phys. Rev. Lett., 109:017202, Jul 2012.
  * [54] Maksym Serbyn, Z. Papić, and Dmitry A. Abanin.  Universal slow growth of entanglement in interacting strongly disordered systems.  Phys. Rev. Lett., 110:260601, Jun 2013.
  * [55] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol.  From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics.  Advances in Physics, 65(3):239–362, 2016.
  * [56] Wilhelm Magnus.  On the exponential solution of differential equations for a linear operator.  Communications on Pure and Applied Mathematics, 7(4):649–673, 1954\.
  * [57] M. Bukov, L. D’Alessio, and A. Polkovnikov.  Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering.  Advances in Physics, 64(2):139–226, 2015.
  * [58] S. Blanes, F. Casas, J.A. Oteo, and J. Ros.  The magnus expansion and some of its applications.  Physics Reports, 470(5):151 – 238, 2009.
  * [59] J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. Daley.  Entanglement growth in quench dynamics with variable range interactions.  Phys. Rev. X, 3:031015, Sep 2013.
  * [60] Maksym Serbyn, Z. Papić, and Dmitry A. Abanin.  Local conservation laws and the structure of the many-body localized states.  Phys. Rev. Lett., 111:127201, Sep 2013.
  * [61] David A. Huse, Rahul Nandkishore, and Vadim Oganesyan.  Phenomenology of fully many-body-localized systems.  Phys. Rev. B, 90:174202, Nov 2014.
  * [62] John Z. Imbrie.  On many-body localization for quantum spin chains.  Journal of Statistical Physics, 163(5):998–1048, 2016.
  * [63] Rajibul Islam, Ruichao Ma, Philipp M. Preiss, M. Eric Tai, Alexander Lukin, Matthew Rispoli, and Markus Greiner.  Measuring entanglement entropy in a quantum many-body system.  Nature, 528(7580):77–83, 2015.
  * [64] Fernando Casas.  Sufficient conditions for the convergence of the magnus expansion.  Journal of Physics A: Mathematical and Theoretical, 40(50):15001–15017, nov 2007.
  * [65] Per Christian Moan and Jitse Niesen.  Convergence of the magnus series.  Foundations of Computational Mathematics, 8(3):291–301, 2008.
  * [66] Norbert Schuch, Michael M. Wolf, Frank Verstraete, and J. Ignacio Cirac.  Entropy scaling and simulability by matrix product states.  Phys. Rev. Lett., 100:030504, Jan 2008.