# A model of quantum gravity on a noisy quantum computer

Muhammad Asaduzzaman  ID<EMAIL_ADDRESS>Department of Physics
and Astronomy, The University of Iowa, Iowa City, Iowa IA 52242, USA  Raghav
G. Jha  ID<EMAIL_ADDRESS>Thomas Jefferson National Accelerator
Facility, Newport News, VA 23606, USA  Bharath Sambasivam  ID<EMAIL_ADDRESS>Department of Physics, Syracuse University, Syracuse, NY 13244, USA

###### Abstract

We study the Sachdev-Ye-Kitaev (SYK) model – an important toy model for
quantum gravity on IBM’s superconducting qubit quantum computers. By using a
graph-coloring algorithm to minimize the number of commuting clusters of terms
in the qubitized Hamiltonian, we find the circuit complexity of the time
evolution using the first-order Lie product formula for $N$ Majorana fermions
is $\mathcal{O}(N^{5}J^{2}t^{2}/\epsilon)$ where $J$ is the dimensionful
coupling parameter, $t$ is the evolution time, and $\epsilon$ is the desired
accuracy. This complexity is a significant improvement over existing result in
the literature. With this improved resource requirement, we perform the time
evolution for $N=6,8$ using up to 300 two-qubit gates and perform different
error mitigation schemes on the noisy hardware results. We find good agreement
with the results obtained using exact diagonalization on classical computers
and noiseless simulators. In particular, we compute the return probability to
the vacuum state after time $t$ and out-of-time order correlators (OTOC) which
is a standard method of quantifying the chaotic nature of quantum many-body
systems.

##  I Introduction

The holographic duality Witten:1998qj  relates a special class of quantum
field theories in $d$ dimensions and quantum gravity in $d+1$ dimensions. This
strong/weak duality enables one to study the properties of strongly coupled
field theory using classical supergravity computations and vice versa.
However, there are no cases where both sides of the duality can be studied
simultaneously. This is mostly due to the lack of analytical methods to deal
with the non-perturbative limit of field theories when they are not
integrable. Several attempts have been made on the lattice in Euclidean time
using Monte Carlo Catterall:2017lub ; Catterall:2020nmn  to study these
theories but they have their own limitations. Therefore, it is often of
interest to find simpler models that have holographic properties and can be
studied in the strong coupling limit. One such model which has now been
extensively studied is the SYK model Sachdev:1992fk ; Kitaev:2015 ;
Sachdev:2015efa ; Maldacena:2016hyu . It is a model of $N$ interacting
Majorana fermions in $0+1$ dimensions with random couplings between $q$
fermions at a time chosen from a Gaussian distribution with zero mean and
variance proportional to $J^{2}/N^{q}$.

An interesting feature of the SYK model is that it develops an approximate
conformal symmetry in the large $N$, low-temperature limit i.e., $N\gg\beta
J\gg 1$ ($\beta$ is the inverse temperature), where it is related to near
extremal black holes that develop the nAdS2 (near AdS2) geometry. It was shown
that this model saturates the chaos bound Maldacena:2015waa , a feature that
is the smoking gun of the holographic (gravity dual) behavior. Since this is a
model with all-to-all fermion interactions in 0+1 dimensions, it is also
computationally easier to handle than other holographic field theories in
higher dimensions. Several groups have extensively studied this model based on
exact diagonalization methods reaching up to 60 Majorana fermions, see Refs.
Gur-Ari:2018okm ; Kobrin:2020xms  and references therein for details. In
addition to the thermodynamic properties, the time evolution complexity in
this model was also studied in Ref. Yang:2019iav  by relating it to the
eigenvalues of the diagonalized Hamiltonian.

As a potential toy model for understanding quantum gravity, it is, therefore,
crucial to study the real-time dynamics of this model beyond methods
accessible by classical computing methods. This leads to the possibility of
studying the SYK model using quantum computers. This direction has already
been explored starting with Ref. Garcia-Alvarez:2016wem . However, they did
not perform the computation on real noisy quantum hardware. In another
interesting work Luo:2017bno , the authors performed time evolution of a
generalized SYK model using a four-qubit NMR quantum simulator and computed
bosonic correlation functions, and studied the late-time behavior. In addition
to the Trotter approach we use in our study, there have been other proposals
for time evolution to improve the circuit complexity (beyond the NISQ era) and
we refer the reader to Ref. Babbush:2018mlj  for details.

In this letter, we study the SYK model on noisy superconducting quantum
computers for the _first_ time and put forth an improved circuit
complexity111Our definition of circuit complexity is the standard one. The
complexity $\mathcal{C}(U)$ is the least number of two-qubit gates in the
circuit that implements the unitary time evolution $U$ of a particular
Hamiltonian $H$. to simulate the SYK model on quantum hardware compared to the
requirements of Ref. Garcia-Alvarez:2016wem . Specifically, we find a
reduction in resource requirement from
$\mathcal{O}(N^{10}J^{2}t^{2}/\epsilon)$ to
$\mathcal{O}(N^{5}J^{2}t^{2}/\epsilon)$ for the Lie-Trotter-based algorithm
2019arXiv191208854C . Using this improvement, we study the time evolution of
this model up to _eight_ Trotter steps on real quantum hardware available
through IBM and compute the return probability and four-point out-of-time
ordered correlators.

##  II SYK Hamiltonian

The Hamiltonian for the SYK model with $N$-Majorana fermions and $q$-fermion
interaction terms is given by:

$H=\frac{(i)^{q/2}}{q!}\sum_{i,j,k,\cdots,q=1}^{N}J_{ijk\cdots
q}~{}~{}\chi_{i}\chi_{j}\chi_{k}\cdots\chi_{q},$ (1)

where $\chi$ are the Majorana fermions satisfying
$\\{\chi_{i},\chi_{j}\\}=\chi_{i}\chi_{j}+\chi_{j}\chi_{i}=\delta_{ij}$. We
consider the $q=4$ case with random all-to-all quartic interactions averaged
over disorder. The random (real) couplings $J_{ijkl}$ are sampled from a
Gaussian distribution with the mean $\overline{J_{ijkl}}=0$ and variance equal
to $\overline{J_{ijkl}^{2}}=\frac{3!J^{2}}{N^{3}}$. We set $J=1$ in this work.
The dimension of the Hilbert space is $\text{dim}(\mathcal{H})=2^{N/2}$ where
$n=N/2$ is the number of qubits222We have $N=2n$ since two Majorana fermions
can be represented by one complex spinless fermion which can be represented by
single qubit. The SYK model can also be considered for $q>4$ interactions and
is exactly solvable in the large $q$ limit Maldacena:2016hyu .

###  II.1 Qubitization and Trotterization

In order to perform the time evolution of the initial state using the SYK
Hamiltonian, we first have to map the fermionic Hamiltonian to qubits. There
are different approaches to this encoding in one and higher dimensions. We
will use the standard method based on the Jordan-Wigner (JW) transformation.
The $N$ fermions, $\chi$ (each of size $2^{n}\times 2^{n})$ in (1) can be
written in terms of tensor product of $N$/2 Pauli matrices $X,Y,Z$ and the
identity matrix $\mathbb{1}$ Garcia-Alvarez:2016wem ; Luo:2017bno  using the
JW map as follows:

$\displaystyle\chi_{2k-1}$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}\prod_{j=1}^{k-1}Z_{j}\Big{)}X_{n}\mathbb{1}^{\otimes(N-2k)/2},$
$\displaystyle\chi_{2k}$
$\displaystyle=\frac{1}{\sqrt{2}}\Big{(}\prod_{j=1}^{k-1}Z_{j}\Big{)}Y_{n}\mathbb{1}^{\otimes(N-2k)/2},$
(2)

where we have suppressed the tensor product between the Pauli matrices to
simplify the notation and the factor in front is to ensure the normalization
of Majorana as discussed in Ref. Kitaev:2015 . In order to simulate the
dynamics on quantum hardware, we first break the Hamiltonian into Pauli
strings as $H=\sum_{j=1}^{m}H_{j}$ and then use the standard Lie-Trotter
formula Seth1996 :

$e^{-iHt}=\Big{(}\prod_{j=1}^{m}e^{-iH_{j}t/r}\Big{)}^{r}+\mathcal{O}\Big{(}\sum_{j<k}\bigg{|}\bigg{|}[H_{j},H_{k}]\bigg{|}\bigg{|}\frac{t^{2}}{r}\Big{)},$
(3)

where we denote the spectral norm by $\left\lVert\cdot\right\rVert$. If the
terms in the decomposition of $H$ are near to commuting, then the Trotter
error is greatly reduced and vanishes if they commute. This is usually not the
case for interesting physical systems. Though the number of terms into which
$H$ is split is $m={N\choose 4}$, one can reduce this number by only summing
over some small number of clusters of Pauli strings $\mathcal{N}\ll m$. This
helps in reducing the Trotter error since the number of terms summed in (3) is
reduced. We find that for $N=6$, the reduction factor i.e., $m/\mathcal{N}=3$
while for $N=8$, $m/\mathcal{N}=35/3$ and this enables us to reliably evolve
to much larger times by controlling the Trotter error considerably (see
Supplemental Material for details). If we demand that the total error in
simulating the time evolution is $\epsilon$, then we need to perform
$r=\mathcal{O}(t^{2}/\epsilon)$ Trotter steps assuming the spectral norm of
commutators to be upper bounded by unity. Using these arguments, Ref. Garcia-
Alvarez:2016wem  estimated the circuit complexity of
$\mathcal{O}(N^{10}t^{2}/\epsilon)$ but this not optimal. The method based on
Lie-Trotter expansion is not the only way of simulating sparse Hamiltonians. A
different way is to use the linear combination of unitaries (LCU) method and
this was done in Ref. Babbush:2018mlj  to bring the cost down to
$\mathcal{O}(N^{7/2}t+N^{5/2}t~{}\text{polylog}(N/\epsilon))$. However, this
is not amenable to NISQ-era hardware implementation. There are also
alternatives to Trotter-based simulations in the NISQ-era Preskill:2018jim
that are more efficient to implement on quantum hardware 2021Quant ; 2022ScPP
but we do not pursue them in this work. We believe it would be interesting to
apply these improved techniques to study $N\geq 8$ fermions.

Let us consider the simplest case of $N=4$. In this case, there is just one
Pauli string and we can write the 2-qubit Hamiltonian as $H=-J_{1234}ZZ$ and
the time evolution circuit Kay:2018huf  is given by Fig. 1.

Figure 1: Circuit implementing one Trotter step of SYK Hamiltonian with $N=4$ Majorana fermions. The definition of the gate is standard: $R_{z}(\theta)=\exp(-i\frac{\theta}{2}Z)$. $N$ | Pauli strings | Clusters | Two-qubit gates
---|---|---|---
4 | 1 | 1 | 2
6 | 15 | 5 | 30
8 | 70 | 6 | 110
10 | 210 | 23 | 498
12 | 495 | 57 | 1504
14 | 1001 | 92 | 3560
16 | 1820 | 116 | 6812
18 | 3060 | 175 | 11962
20 | 4845 | 246 | 19984
Table 1: The two-qubit gate cost for the time evolution of the SYK Hamiltonian
for one Trotter step, number of Pauli strings, and the number of clusters
$\mathcal{N}$ of commuting Pauli strings for different $N$.

The two-qubit cost of simulating various $N$ and the number of terms in the
decomposition of $H$ can be found in Table 1 for $N\leq 20$. The circuit
complexity estimate based on commuting Pauli strings exploiting the graph-
coloring algorithm grows333For small $N$ i.e., $N\leq 10$, the two-qubit cost
based on quantum Shannon decomposition (QSD) of the unitary matrix
$\exp(-iHt)$ is lower Shende2006 . However, it is known to have the
exponential scaling for larger $N$. like $\mathcal{O}(N^{5})$, a substantial
improvement over $\mathcal{O}(N^{10})$ proposed in Ref. Garcia-Alvarez:2016wem
. To obtain this improved scaling, we used well-known methods such as
partitioning of the Pauli terms in commuting clusters (using the degree of
saturation at each node i.e., DSatur 10.1145/359094.359101  graph-coloring
algorithm) and then constructing the diagonalizing operator (circuit) for each
cluster van_den_Berg_2020 ; 2020arXiv200105983G ; Miller:2022sol ;
kurita2022pauli ; Murairi:2023oti  separately followed by exponentiation for
time evolution. At the end we add the gate costs from all clusters Jena2019 ;
Verteletskyi_2020 ; Murairi:2022zdg . The number of Pauli strings in the
decomposition of $H$ is given by ${N\choose
4}=N^{4}/24-N^{3}/4+11N^{2}/24-N/4$ and grows like $\mathcal{O}(N^{4}/4!)$ in
the large $N$ limit, and the cost to simulate the simplest non-trivial Pauli
string of length $N/2$ is $\mathcal{O}(N)$. Hence, it appears that the $\sim
N^{5}$ circuit complexity is close to _optimal_ for this approach to
Hamiltonian simulation444For $q$-local SYK Hamiltonian, the complexity will be
$\sim\mathcal{O}(N^{q+1})$ (up to factor of $q!$) and this suggests that
beyond some critical $q$, there won’t be any quantum advantage for simulation
of this model.. The complexity can be improved by using Bravyi-Kitaev mapping
from fermions, however, such optimizations are not required for small $N$
which are accessible in the NISQ era. We provide details of how we estimate
the gate resources in the Supplemental Material (SM) for interested readers.

###  II.2 Return probability

One of the main motivations for using quantum computers is to understand the
real-time evolution of complicated quantum many-body systems. An observable we
compute in this work is the return probability by evolving an initial
state555The initial state chosen here belongs to the set of common eigenstates
of the SYK spin operators defined in Numasawa:2019gnl  and forms a complete
basis of the SYK Hilbert space given by $\ket{\psi_{0}}=\ket{0}^{\otimes n}$
with the SYK Hamiltonian for $N_{t}$ Trotter steps and computing the overlap
as:

$P_{0}=|\langle\psi_{0}|e^{-iHt}|\psi_{0}\rangle|^{2}.$ (4)

The return probability can equivalently be thought of as computing the
expectation value of the projection operator $P_{0}$ onto the $\ket{\psi_{0}}$
state. The return probability is closely related (proportional up to $|t|$ in
the semi-classical regime) to the spectral form factor Liu:2016rdi  and shows
similar behavior of slope, dip, ramp, and plateau upon averaging over an
ensemble of Hamiltonians Cotler:2016fpe ; Numasawa:2019gnl  which are key
features of the SYK model.

The error-mitigated hardware results of the return probabilities for five
realizations of the SYK Hamiltonian with $N=6$ and disorder average are shown
in Fig. 2. The dashed black lines are the results from the exact time
evolution, while the dashed blue curve is the ensemble average of the exact
evolution over the five realizations of the model. As $N$ increases, the
resource requirements quickly increase as shown in Table 1. Therefore, for
$N=8$ we consider the time-evolution of just a single instance of the SYK
Hamiltonian in Fig. 3. The markers are results obtained from the 127-qubit IBM
hardware with various degrees of error-mitigation applied.

Figure 2: The return probability computed for five realisations (panels $a$
through $e$) of the SYK model with $N=6$. The error bars are smaller than the
symbol size. The disordered average (panel $f$) clearly shows the slope region
and the starting of the ramp behavior. The hardware results for the last
Trotter step ($t=12)$ were obtained on ibm_cusco while all others were
obtained on ibm_nazca. For each data point on the IBM device, we ran 75
physics circuits and 75 self-mitigation circuits with 2048 shots. Figure 3:
The return prob computed for a single instance of the SYK model with $N=8$.
The hardware results were obtained on ibm_kyoto. The obtained results for the
return probability are consistent with $1-p$ where $p$ corresponds to the
value when the four-qubit channel (gray dashed line) is completely
depolarized.

We used devices such as ibm_cusco, ibm_nazca and ibm_kyoto. These machines use
the eagle r3 processor, where the native two-qubit gate is not the standard
and familiar CX but the echoed cross-resonance (ECR) gate (see (17) for the
definition of the gate). ECR gates are more resilient against error drifting,
thereby reducing the need for frequent calibration. The leading source of gate
noise in current devices is the two-qubit gates. To deal with this, the error
mitigation strategy we employ is a combination of Pauli twirling/randomized
compiling Wallman_2016  for the ECR gates and self-mitigation ARahman:2022tkr
. For the ECR gate, there is a set of 16 single-qubit conjugations (see SM for
details) that leave the exact implementation of the ECR gate invariant, up to
a global phase. For each gate implementation, the noise is ‘twirled’
differently for these 16 conjugations. This helps convert a majority of the
coherent noise in the circuit into incoherent noise that is well-modeled by
the depolarizing quantum channel (see SM for details).

We then construct mitigation circuits with the same structure and gate count
as the physics circuits, to characterize the probability of errors $p$ due to
the depolarizing channel. With self-mitigation, the idea is to construct a
circuit consisting of $N_{t}/2$ Trotter steps forward in time and $N_{t}/2$
Trotter steps backward in time. On a noiseless device, this would bring one
back to the initial $\ket{\psi_{0}}$ state, up to Trotter error. However, on a
noisy device, the state of the system after the mitigation circuit would be
some $\overline{\ket{\psi_{0}}}$. The probability of error is then
$p=|\bra{\psi_{0}}\overline{{\psi_{0}}\rangle}|^{2}$. Using this, we can
extract the noiseless return probability $\langle P_{0}\rangle$ from the noisy
one $\overline{\langle P_{0}\rangle}$ using PhysRevLett.127.270502

$\langle P_{0}\rangle=\frac{\overline{\langle P_{0}\rangle}-{2^{-n}}p}{1-p}.$
(5)

We build 75 Pauli-twirled physics circuits, and 75 Pauli-twirled self-
mitigation circuits and run them each for 2048 shots. We also use the standard
M3 protocol PRXQuantum.2.040326  for qubit measurement errors and dynamical
decoupling PhysRevA.58.2733 ; ZANARDI199977 ; PhysRevA.59.4178 ;
Ezzell:2022uat  to suppress decoherence noise. The results for the return
probability666We only show even Trotter steps because for self-mitigation, one
requires forward and backward Trotter evolution by an equal amount with
different types of error mitigation applied and exact evolution are shown in
panels (_a-e_) of Fig. 2 for five different instances of the SYK Hamiltonian
with $N=6$ and the disorder average over these realizations is shown in the
last panel (_f_). The return probability for $N=8$ is shown in Fig. 3 for a
single instance on ibm_kyoto. The obtained results are consistent with pure
noise and therefore we did not study different realizations of the model for
this case.

For $N=6$, the self-mitigated return probability matches the results from
exact diagonalization very well even for a relatively large evolution time,
and correspondingly large number of ECR gates in the circuit. For this case,
we are able to probe the slope and dip regions robustly and are beginning to
see signs of the ramp region. To simulate the SYK model with $N=6$ up to
$t=12$ (eight Trotter steps), we require about 300 ECR gates, which is close
to the limit of the current state-of-the-art. For $N=8$, we need about 340 ECR
gates for doing two Trotter steps given the restrictions from qubit topology
on IBM hardware. It appears to be beyond the current capability even with
advanced error mitigation methods.

##  III OTOC computation

An important feature of the SYK model is that for large $N$ and in the low-
temperature limit, it is maximally chaotic and is a quintessential example of
a fast scrambler. A defining feature of such systems is that the quantum
information shared between a small number of elementary degrees of freedom is
rapidly distributed into exponentially many degrees of freedom. This is known
as ‘scrambling’. Since black holes are known to be the fastest scramblers in
nature, finding similar behavior in a non-gravitational quantum many-body
system is highly non-trivial. This is one reason why the SYK model is
considered to be important for understanding aspects of quantum gravity. To
quantify this chaos in a quantum system, one considers the out-of-time (OTO)
commutator between two operators $W$ and $V$ given by Hosur:2015ylk ;
Fu:2016yrv ; Kitaev:2017awl ; Lantagne-Hurtubise:2019svg :

$C(t)=-\langle[W(t),V(0)]^{\dagger}\,[W(t),V(0)]\rangle,$ (6)

In general $C(t)$ starts from zero, and becomes significant at some later time
$t$, which one refers to as scrambling time (also Ehrenfest time) denoted by
$t_{*}$. By considering the commutator expansion, we can define the out-of-
time order correlators (OTOC)777These were first introduced with reference to
the study of disordered superconductors 1969JETP…28.1200L . In fact, OTOC is
also closely related to the thermal average of signals from Loschmidt echo as
explored in PhysRevLett.124.160603  as:

$\text{OTOC}:=F(t)=\langle W(t)V(0)W(t)V(0)\rangle_{\beta},$ (7)

where $W$ and $V$ are generic Hermitian operators and do not have the same
symmetry as the Hamiltonian. The growth of the $C(t)$ is related to the decay
of the OTOCs i.e., $F(t)$ through the simple relation, $C(t)=2(1-F(t))$. The
time-evolution of the operator $W$ with Hamiltonian $H$ in the Heisenberg
representation is $W(t)=W_{t}=e^{iHt}W(0)e^{-iHt}$ and
$\langle\cdot\rangle_{\beta}=\mathrm{Tr}\\{\rho~{}\cdot\\}$ denotes the
thermal average with the density matrix $\rho$ at inverse temperature $\beta$.
A common choice for $\rho$ is to just use the $T\to\infty$ limit given by the
normalized identity matrix, $\mathbb{1}/\text{dim}(\mathcal{H})$
Sunderhauf:2019djv  such that $\rm{Tr}(\rho)=1$. For the $W$ and $V$ we can
either take Pauli matrices or Majorana fermions such that $W=\chi_{i}$ and
$V=\chi_{j}$ with $i\neq j$ and average over the different pairs $(i,j)$. We
used the simplest case of a single Pauli matrix, i.e., $W=V=Z$. In what
follows, we denote $1-F(t)$ by $1-F(t)_{ij}=O_{ij}(t)$ where $i$ and $j$
denote the qubit location of the single-qubit operator $W$ and $V$
respectively and compute this on the hardware. The choice of these operators
does not change the basic features of exponential growth and saturation.

The computation of OTOC in quantum many-body systems is useful because it
allows for the basis of comparison between black holes and quantum many-body
systems without any additional requirements on how the dual gravity theory is
fully defined. The OTOC grows as $\exp(\lambda\tilde{t})$ after all two-point
correlations (i.e., $\langle W(t)W(0)\rangle$) have decayed (often denoted as
$t_{d}$) i.e., the exponential growth occurs for $t_{d}<\tilde{t}\leq t_{*}$.
The exponent $\lambda$ is known as the Lyapunov exponent. As $t$ increases
beyond $t_{*}$, the exponential growth of OTOC saturates. The scrambling
behavior is often a defining feature of quantum many-body systems with all-to-
all (non-local) interactions. In Ref. Maldacena:2015waa , the authors showed
that $\lambda\leq 2\pi/\beta$ for any quantum system. The SYK model exactly
saturates this bound at low temperatures. These observables have been
extensively studied using classical computing methods such as exact
diagonalization Cao:2020rhe ; Anegawa:2023vxq . The current resources do not
allow us to access very large values of $N$ or explore finite values of
$\beta$, but we take a first step at computing these quantities in the
simplest setting on the quantum computer in this work. With the available
resources, we have been able to simulate four Trotter steps getting close to
the scrambling time888All time scales are in units where $J=1$ (which we
estimate to be around $t\sim 10$ for $N=6$) on quantum hardware, with an
improved hardware resource requirement.

Figure 4: The OTOC for three realisations (panels $a$–$c$) of the SYK model
with $N=6$. The disordered average is shown in panel ($d$). Hardware results
are obtained on ibm_kyoto, and ibm_cusco depending on the availability of the
device during the run time. $N_{u}=600$ unitaries were used upto first Trotter
step and at most $N_{u}=900$ unitaries were used for the rest of the
measurements of $\langle W\rangle$ and $\langle VW(t)V\rangle$ with
$N_{M}=4000$ shots per circuit.

###  III.1 Results for OTOC

Even though the OTOC seems simple to compute, the experimental/quantum
computer measurement of the OTOC is very challenging because of the unusual
time ordering. In the first study of its kind, Ref. PhysRevX.7.031011  studied
OTOC of local operators for the first time on a nuclear magnetic resonance
(NMR) quantum simulator, and soon after it was measured on trapped-ion quantum
magnet Garttner:2016mqj . In addition to these efforts, to compute OTOC on
quantum hardware, we need to define a protocol for performing this
measurement. Several proposals have been put forth Yao:2016ayk ;
2018arXiv180600781A ; Vermersch:2018sru ; Mi:2021gdf  in this regard.
Recently, the protocol proposed by Vermersch et al. has been used to compute
modified OTOCs in the Ising spin model on hyperbolic space Asaduzzaman:2023htk
on IBM hardware. We discuss the details of the global protocol (the one we
use) based on randomized measurements proposed in Ref. Vermersch:2018sru  in
the SM. The protocol measures two independent quantities and from the
correlation of the measurements, normalized out of time ordered correlators
are computed using the following definition

$O_{ij}(t)=\frac{\overline{\langle W_{i}(t)\rangle_{u,k_{0}}\langle
V_{j}^{\dagger}W_{i}(t)V_{j}\rangle_{u,k_{0}}}}{\overline{\langle
W_{i}(t)\rangle_{u,k_{0}}^{2}}}.$ (8)

The $\overline{\cdots}$ denotes the ensemble average of an operator over a set
of unitaries $\\{u,v,\cdots\\}$ and $\langle\cdots\rangle_{u,k_{0}}$ denotes
quantum expectation values computed from measurements with a particular
unitary $u$ and a chosen initial state $\ket{k_{0}}$.

The protocol introduces computational overhead due to running circuits for
many different unitaries $N_{u}$ for the computation of each data point.
However, there are some salient features of this protocol that make this a
suitable choice. The normalized quantity in (8) consists of a ratio, and hence
the protocol is robust against depolarizing error and readout errors since
this scales the estimated quantities $\langle W(t)\rangle$ and $\langle
VW(t)V\rangle$ at the same rate. This resembles the Twirled readout error
extinction (T-Rex) error mitigation technique van2022model  which involves
twirling the circuit with Pauli matrices at the end before measurement. The
protocol is also robust against imperfect creation of the random unitaries
Vermersch:2018sru . It remains to be seen if the protocol can be made more
cost-effective by designing the unitaries with approximate construction like
that of unitary $t$-design dankert2009exact ; Nakata:2021tmt . The protocol
does not require any ancilla qubits and the backward propagation is not
required making it viable for the implementation with the NISQ devices. An
extensive study comparing different protocols for the SYK model on hardware is
an interesting direction on which we hope to report in the future. The results
for OTOC from ibm_cusco and ibm_kyoto for three instances of $N=6$ and the
disorder average over them is shown in Fig. 4. For this computation, we used
the M3 readout error mitigation protocol and dynamical decoupling. Even
without the self mitigation, we find good agreement with the exact results. We
believe this is due to the robust nature of the protocol against different
error channels. We conclude the section, by commenting that there will be
statistical errors associated with finite $N_{u}$ and $N_{M}$. In principle,
with infinite measurements $N_{M}$, the overall error is saturated to
$\frac{1}{\sqrt{N_{u}}}$ in noiseless environment Vermersch:2018sru . Since
with finite resources, this is not viable, it is imperative to run circuits
with different number of shots $N_{M}$ and observe the error rate by gradually
increasing the number of unitaries $N_{u}$.

##  IV Summary and Discussion

We have proposed an improved circuit complexity for Hamiltonian simulation of
the SYK model with $N$ Majorana fermions from
$\mathcal{O}(N^{10}t^{2}/\epsilon)$ down to $\mathcal{O}(N^{5}t^{2}/\epsilon)$
using well-known techniques and performed quantum simulations of the SYK model
on noisy quantum computers. We studied the behavior of return probability due
to its close relation to the spectral form factor in the semi-classical limit
using 127-qubit ibm quantum computers for $N=6,8$ Majorana fermions. Due to
the noisy devices currently accessible for the computations, we applied state-
of-the-art error mitigation methods to the hardware result and showed that it
agrees well with the exact time evolution results.

We also computed the out-of-time order correlators which are diagnostic of the
chaotic behavior of quantum many-body systems for $N=6$. Though we cannot
study the large $N$ and large $\beta$ limits and see signs of saturation of
the Lyapunov exponent on current NISQ devices, we believe that our work will
be useful in future explorations of this model beyond the NISQ era. In this
regard, to study quantum chaos on a quantum computer, it might also be useful
to study OTOC in simplified models similar to SYK Xu:2020shn ; Garcia-
Garcia:2020cdo ; Hanada:2023rkf  that are conjectured to have the same low-
temperature, large $N$ dynamics as the pure SYK model we considered here.
Another direction for extension of our work is to consider $q>4$ and explore
the resource requirements and time evolution in detail. It would be very
useful to study the dynamics of the model over different time scales extending
to the plateau region for the return probability at finite $\beta$. These
interesting problems would require resources that are beyond the contemporary
NISQ-era hardware. We therefore leave these questions for future work. The use
of quantum computers for models such as the SYK model in coming decades will
not only provide new insights into the holographic principle but also into the
interesting world of strange metals and quantum many-body systems
Chowdhury:2021qpy .

### Acknowledgements

MA is supported under the U.S. Department of Energy grant DE-SC0019139. RGJ is
supported by the U.S. Department of Energy, Office of Science, National
Quantum Information Science Research Centers, Co-design Center for Quantum
Advantage (C2QA) under contract number DE-SC0012704 and by the U.S. Department
of Energy, Office of Science, Office of Nuclear Physics under contract number
DE-AC05-06OR23177. BS is supported in part by the U.S. Department of Energy,
Office of Science, Office of High Energy Physics, under Award Number DE-
SC0009998. MA and BS would like to thank Jefferson Lab’s Quantum Computing
Bootcamp for the hospitality where this work started. We thank the IBM quantum
hub at Brookhaven National Laboratory for providing access to the IBMQ quantum
computers on which the computations were done.

###  IV.1 Data Availability Statement

The $N=6,8$ SYK Hamiltonian realizations, its decomposition in terms of Pauli
strings, and the time evolution circuit as open qasm 2.0 files for single
Trotter step can be obtained from Ref. Asad_SYK_2023 .

## References

  * (1) E. Witten, “Anti-de Sitter space and holography,” Adv. Theor. Math. Phys. 2 (1998) 253–291, arXiv:hep-th/9802150.
  * (2) S. Catterall, R. G. Jha, D. Schaich, and T. Wiseman, “Testing holography using lattice super-Yang-Mills theory on a 2-torus,” Phys. Rev. D 97 no. 8, (2018) 086020, arXiv:1709.07025 [hep-th].
  * (3) S. Catterall, J. Giedt, R. G. Jha, D. Schaich, and T. Wiseman, “Three-dimensional super-Yang–Mills theory on the lattice and dual black branes,” Phys. Rev. D 102 no. 10, (2020) 106009, arXiv:2010.00026 [hep-th].
  * (4) S. Sachdev and J. Ye, “Gapless spin fluid ground state in a random, quantum Heisenberg magnet,” Phys. Rev. Lett. 70 (1993) 3339, arXiv:cond-mat/9212030.
  * (5) A. Kitaev, “A simple model of quantum holography,”. http://online.kitp.ucsb.edu/online/entangled15/.
  * (6) S. Sachdev, “Bekenstein-Hawking Entropy and Strange Metals,” Phys. Rev. X 5 no. 4, (2015) 041025, arXiv:1506.05111 [hep-th].
  * (7) J. Maldacena and D. Stanford, “Remarks on the Sachdev-Ye-Kitaev model,” Phys. Rev. D 94 no. 10, (2016) 106002, arXiv:1604.07818 [hep-th].
  * (8) J. Maldacena, S. H. Shenker, and D. Stanford, “A bound on chaos,” JHEP 08 (2016) 106, arXiv:1503.01409 [hep-th].
  * (9) G. Gur-Ari, R. Mahajan, and A. Vaezi, “Does the SYK model have a spin glass phase?,” JHEP 11 (2018) 070, arXiv:1806.10145 [hep-th].
  * (10) B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, “Many-Body Chaos in the Sachdev-Ye-Kitaev Model,” Phys. Rev. Lett. 126 no. 3, (2021) 030602, arXiv:2002.05725 [hep-th].
  * (11) R.-Q. Yang and K.-Y. Kim, “Time evolution of the complexity in chaotic systems: a concrete example,” JHEP 05 (2020) 045, arXiv:1906.02052 [hep-th].
  * (12) L. García-Álvarez, I. L. Egusquiza, L. Lamata, A. del Campo, J. Sonner, and E. Solano, “Digital Quantum Simulation of Minimal AdS/CFT,” Phys. Rev. Lett. 119 no. 4, (2017) 040501, arXiv:1607.08560 [quant-ph].
  * (13) Z. Luo, Y.-Z. You, J. Li, C.-M. Jian, D. Lu, C. Xu, B. Zeng, and R. Laflamme, “Quantum Simulation of the Non-Fermi-Liquid State of Sachdev-Ye-Kitaev Model,” npj Quantum Inf. 5 (2019) 53, arXiv:1712.06458 [quant-ph].
  * (14) R. Babbush, D. W. Berry, and H. Neven, “Quantum Simulation of the Sachdev-Ye-Kitaev Model by Asymmetric Qubitization,” Phys. Rev. A 99 no. 4, (2019) 040301, arXiv:1806.02793 [quant-ph].
  * (15) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, “A Theory of Trotter Error,” arXiv e-prints (Dec., 2019) arXiv:1912.08854, arXiv:1912.08854 [quant-ph].
  * (16) S. Lloyd, “Universal quantum simulators,” Science 273 no. 5278, (Aug., 1996) 1073–1078. https://doi.org/10.1126/science.273.5278.1073.
  * (17) J. Preskill, “Quantum Computing in the NISQ era and beyond,” Quantum 2 (2018) 79, arXiv:1801.00862 [quant-ph].
  * (18) S. Barison, F. Vicentini, and G. Carleo, “An efficient quantum algorithm for the time evolution of parameterized circuits,” Quantum 5 (July, 2021) 512, arXiv:2101.04579 [quant-ph].
  * (19) J. W. Z. Lau, T. Haug, L. C. Kwek, and K. Bharti, “NISQ Algorithm for Hamiltonian simulation via truncated Taylor series,” SciPost Physics 12 no. 4, (Apr., 2022) 122, arXiv:2103.05500 [quant-ph].
  * (20) A. Kay, “Tutorial on the Quantikz Package,” arXiv:1809.03842 [quant-ph].
  * (21) V. Shende, S. Bullock, and I. Markov, “Synthesis of quantum-logic circuits,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25 no. 6, (June, 2006) 1000–1010. https://doi.org/10.1109/tcad.2005.855930.
  * (22) D. Brélaz, “New methods to color the vertices of a graph,” Commun. ACM 22 no. 4, (Apr, 1979) 251–256. https://doi.org/10.1145/359094.359101.
  * (23) E. van den Berg and K. Temme, “Circuit optimization of hamiltonian simulation by simultaneous diagonalization of pauli clusters,” Quantum 4 (Sep, 2020) 322. https://doi.org/10.22331%2Fq-2020-09-12-322.
  * (24) K. Gui, T. Tomesh, P. Gokhale, Y. Shi, F. T. Chong, M. Martonosi, and M. Suchara, “Term Grouping and Travelling Salesperson for Digital Quantum Simulation,” arXiv e-prints (Jan., 2020) arXiv:2001.05983, arXiv:2001.05983 [quant-ph].
  * (25) D. Miller, L. E. Fischer, I. O. Sokolov, P. K. Barkoutsos, and I. Tavernelli, “Hardware-Tailored Diagonalization Circuits,” arXiv:2203.03646 [quant-ph].
  * (26) T. Kurita, M. Morita, H. Oshima, and S. Sato, “Pauli string partitioning algorithm with the ising model for simultaneous measurement,” 2022.
  * (27) E. M. Murairi and M. J. Cervia, “Reducing Circuit Depth with Qubitwise Diagonalization,” arXiv:2306.00170 [quant-ph].
  * (28) A. Jena, S. Genin, and M. Mosca, “Pauli Partitioning with Respect to Gate Sets,” arXiv e-prints (July, 2019) arXiv:1907.07859, arXiv:1907.07859 [quant-ph].
  * (29) V. Verteletskyi, T.-C. Yen, and A. F. Izmaylov, “Measurement optimization in the variational quantum eigensolver using a minimum clique cover,” J. Chem. Phys. 152 no. 12, (Mar., 2020) 124114, arXiv:1907.03358 [quant-ph].
  * (30) E. M. Murairi, M. J. Cervia, H. Kumar, P. F. Bedaque, and A. Alexandru, “How many quantum gates do gauge theories require?,” Phys. Rev. D 106 no. 9, (2022) 094504, arXiv:2208.11789 [hep-lat].
  * (31) T. Numasawa, “Late time quantum chaos of pure states in random matrices and in the Sachdev-Ye-Kitaev model,” Phys. Rev. D 100 no. 12, (2019) 126017, arXiv:1901.02025 [hep-th].
  * (32) Y. Liu, M. A. Nowak, and I. Zahed, “Disorder in the Sachdev-Yee-Kitaev Model,” Phys. Lett. B 773 (2017) 647–653, arXiv:1612.05233 [hep-th].
  * (33) J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher, and M. Tezuka, “Black Holes and Random Matrices,” JHEP 05 (2017) 118, arXiv:1611.04650 [hep-th]. [Erratum: JHEP 09, 002 (2018)].
  * (34) J. J. Wallman and J. Emerson, “Noise tailoring for scalable quantum computation via randomized compiling,” Phys. Rev. A 94 no. 5, (Nov., 2016) 052325, arXiv:1512.01098 [quant-ph].
  * (35) S. A Rahman, R. Lewis, E. Mendicelli, and S. Powell, “Self-mitigating Trotter circuits for SU(2) lattice gauge theory on a quantum computer,” Phys. Rev. D 106 no. 7, (2022) 074502, arXiv:2205.09247 [hep-lat].
  * (36) M. Urbanek, B. Nachman, V. R. Pascuzzi, A. He, C. W. Bauer, and W. A. de Jong, “Mitigating depolarizing noise on quantum computers with noise-estimation circuits,” Phys. Rev. Lett. 127 (Dec, 2021) 270502. https://link.aps.org/doi/10.1103/PhysRevLett.127.270502.
  * (37) P. D. Nation, H. Kang, N. Sundaresan, and J. M. Gambetta, “Scalable mitigation of measurement errors on quantum computers,” PRX Quantum 2 (Nov, 2021) 040326. https://link.aps.org/doi/10.1103/PRXQuantum.2.040326.
  * (38) L. Viola and S. Lloyd, “Dynamical suppression of decoherence in two-state quantum systems,” Phys. Rev. A 58 (Oct, 1998) 2733–2744. https://link.aps.org/doi/10.1103/PhysRevA.58.2733.
  * (39) P. Zanardi, “Symmetrizing evolutions,” Physics Letters A 258 no. 2, (1999) 77–82. https://www.sciencedirect.com/science/article/pii/S0375960199003655.
  * (40) D. Vitali and P. Tombesi, “Using parity kicks for decoherence control,” Phys. Rev. A 59 (Jun, 1999) 4178–4186. https://link.aps.org/doi/10.1103/PhysRevA.59.4178.
  * (41) N. Ezzell, B. Pokharel, L. Tewala, G. Quiroz, and D. A. Lidar, “Dynamical decoupling for superconducting qubits: a performance survey,” arXiv:2207.03670 [quant-ph].
  * (42) P. Hosur, X.-L. Qi, D. A. Roberts, and B. Yoshida, “Chaos in quantum channels,” JHEP 02 (2016) 004, arXiv:1511.04021 [hep-th].
  * (43) W. Fu and S. Sachdev, “Numerical study of fermion and boson models with infinite-range random interactions,” Phys. Rev. B 94 no. 3, (2016) 035135, arXiv:1603.05246 [cond-mat.str-el].
  * (44) A. Kitaev and S. J. Suh, “The soft mode in the Sachdev-Ye-Kitaev model and its gravity dual,” JHEP 05 (2018) 183, arXiv:1711.08467 [hep-th].
  * (45) E. Lantagne-Hurtubise, S. Plugge, O. Can, and M. Franz, “Diagnosing quantum chaos in many-body systems using entanglement as a resource,” Phys. Rev. Res. 2 no. 1, (2020) 013254, arXiv:1907.01628 [cond-mat.str-el].
  * (46) A. I. Larkin and Y. N. Ovchinnikov, “Quasiclassical Method in the Theory of Superconductivity,” Soviet Journal of Experimental and Theoretical Physics 28 (June, 1969) 1200.
  * (47) B. Yan, L. Cincio, and W. H. Zurek, “Information scrambling and loschmidt echo,” Phys. Rev. Lett. 124 (Apr, 2020) 160603. https://link.aps.org/doi/10.1103/PhysRevLett.124.160603.
  * (48) C. Sünderhauf, L. Piroli, X.-L. Qi, N. Schuch, and J. I. Cirac, “Quantum chaos in the Brownian SYK model with large finite $N$: OTOCs and tripartite information,” JHEP 11 (2019) 038, arXiv:1908.00775 [quant-ph].
  * (49) Y. Cao, Y.-N. Zhou, T.-T. Shi, and W. Zhang, “Towards quantum simulation of Sachdev-Ye-Kitaev model,” Sci. Bull. 65 (2020) 1170–1176, arXiv:2003.01514 [cond-mat.dis-nn].
  * (50) T. Anegawa, N. Iizuka, A. Mukherjee, S. K. Sake, and S. P. Trivedi, “Sparse random matrices and Gaussian ensembles with varying randomness,” arXiv:2305.07505 [hep-th].
  * (51) J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, “Measuring out-of-time-order correlators on a nuclear magnetic resonance quantum simulator,” Phys. Rev. X 7 (Jul, 2017) 031011. https://link.aps.org/doi/10.1103/PhysRevX.7.031011.
  * (52) M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, “Measuring out-of-time-order correlations and multiple quantum spectra in a trapped ion quantum magnet,” Nature Phys. 13 (2017) 781, arXiv:1608.08938 [quant-ph].
  * (53) N. Y. Yao, F. Grusdt, B. Swingle, M. D. Lukin, D. M. Stamper-Kurn, J. E. Moore, and E. A. Demler, “Interferometric Approach to Probing Fast Scrambling,” arXiv:1607.01801 [quant-ph].
  * (54) D. Aggarwal, S. Raj, B. K. Behera, and P. K. Panigrahi, “Application of quantum scrambling in Rydberg atom on IBM quantum computer,” arXiv e-prints (June, 2018) arXiv:1806.00781, arXiv:1806.00781 [quant-ph].
  * (55) B. Vermersch, A. Elben, L. M. Sieberer, N. Y. Yao, and P. Zoller, “Probing scrambling using statistical correlations between randomized measurements,” Phys. Rev. X 9 no. 2, (2019) 021061, arXiv:1807.09087 [quant-ph].
  * (56) X. Mi et al., “Information scrambling in quantum circuits,” Science 374 no. 6574, (2021) abg5029, arXiv:2101.08870 [quant-ph].
  * (57) M. Asaduzzaman, S. Catterall, Y. Meurice, and G. C. Toga, “Quantum Ising model on two dimensional anti-de Sitter space,” arXiv:2309.04383 [quant-ph].
  * (58) E. Van Den Berg, Z. K. Minev, and K. Temme, “Model-free readout-error mitigation for quantum expectation values,” Physical Review A 105 no. 3, (2022) 032620.
  * (59) C. Dankert, R. Cleve, J. Emerson, and E. Livine, “Exact and approximate unitary 2-designs and their application to fidelity estimation,” Physical Review A 80 no. 1, (2009) 012304.
  * (60) Y. Nakata et al., “Quantum Circuits for Exact Unitary t-Designs and Applications to Higher-Order Randomized Benchmarking,” PRX Quantum 2 no. 3, (2021) 030339, arXiv:2102.12617 [quant-ph].
  * (61) S. Xu, L. Susskind, Y. Su, and B. Swingle, “A Sparse Model of Quantum Holography,” arXiv:2008.02303 [cond-mat.str-el].
  * (62) A. M. García-García, Y. Jia, D. Rosa, and J. J. M. Verbaarschot, “Sparse Sachdev-Ye-Kitaev model, quantum chaos and gravity duals,” Phys. Rev. D 103 no. 10, (2021) 106002, arXiv:2007.13837 [hep-th].
  * (63) M. Hanada, A. Jevicki, X. Liu, E. Rinaldi, and M. Tezuka, “A model of randomly-coupled Pauli spins,” arXiv:2309.15349 [hep-th].
  * (64) D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, “Sachdev-Ye-Kitaev models and beyond: Window into non-Fermi liquids,” Rev. Mod. Phys. 94 no. 3, (2022) 035004, arXiv:2109.05037 [cond-mat.str-el].
  * (65) M. Asaduzzaman, R. G. Jha, and B. Sambasivam, “A model of quantum gravity on a noisy quantum computer - notebooks and quantum circuits release,” Zenodo (November, 2023) . https://zenodo.org/record/10202045.
  * (66) J. K. Iverson and J. Preskill, “Coherence in logical quantum channels,” New Journal of Physics 22 no. 7, (July, 2020) 073066, arXiv:1912.04319 [quant-ph].
  * (67) M. Urbanek, B. Nachman, V. R. Pascuzzi, A. He, C. W. Bauer, and W. A. de Jong, “Mitigating depolarizing noise on quantum computers with noise-estimation circuits,” Phys. Rev. L 127 (2021) 270502, arXiv:2103.08591 [quant-ph].
  * (68) A. He, B. Nachman, W. A. de Jong, and C. W. Bauer, “Zero-noise extrapolation for quantum-gate error mitigation with identity insertions,” Phys. Rev. A 102 no. 1, (2020) 012426, arXiv:2003.04941 [quant-ph].
  * (69) P. Nation, “Generating Pauli-twirled circuits in Qiskit,”. https://quantum-enablement.org/posts/2023/2023-02-02-pauli_twirling.html.
  * (70) X. Chen, T. Zhou, and C. Xu, “Measuring the distance between quantum many-body wave functions,” Journal of Statistical Mechanics: Theory and Experiment 2018 no. 7, (2018) 073101.
  * (71) M. Asaduzzaman, S. Catterall, G. C. Toga, Y. Meurice, and R. Sakai, “Quantum simulation of the N-flavor Gross-Neveu model,” Phys. Rev. D 106 no. 11, (2022) 114515, arXiv:2208.05906 [hep-lat].
  * (72) D. C. McKay, I. Hincks, E. J. Pritchett, M. Carroll, L. C. G. Govia, and S. T. Merkel, “Benchmarking Quantum Processor Performance at Scale,” arXiv:2311.05933 [quant-ph].

##  Appendix A Supplemental Material (SM) for ‘A model of quantum gravity on a
noisy quantum computer’

## Optimizing the time evolution circuit and return probability

One of the main results of the paper is how to perform efficient time
evolution of the model. We now provide details related to optimizing the gate
costs of the time evolution on the basis of commuting Pauli strings in the
decomposition of the qubit Hamiltonian. Let us take $N=6$ as an example. In
this case, we have a three-qubit Hamiltonian which can be decomposed into 15
Pauli strings each with random coupling taken from a well-defined
distribution. The generic Hamiltonian is written as:

$\displaystyle
H=\alpha_{1}\mathbb{1}ZZ+\alpha_{2}X\mathbb{1}X+\alpha_{3}X\mathbb{1}Y+\alpha_{4}XXZ+\alpha_{5}XYZ+\alpha_{6}Y\mathbb{1}X+\alpha_{7}Y\mathbb{1}Y+\alpha_{8}YXZ+\alpha_{9}YYZ+\alpha_{10}Z\mathbb{1}Z$
(9) $\displaystyle+\alpha_{11}ZXX+\alpha_{12}ZXY\newline
+\alpha_{13}ZYX+\alpha_{14}ZYY+\alpha_{15}ZZ\mathbb{1},$

where $\alpha_{1,2,\dots,15}$ are the random couplings for a single
realization of the model. By inspection, we find that it is possible to divide
the terms of $H$ into five clusters where all Pauli strings in a given cluster
commute with each other. The clusters are:

$\\{\mathbb{1}ZZ,ZXX,ZYY\\},\\{X\mathbb{1}X,YXZ,ZXY\\},\\{X\mathbb{1}Y,Y\mathbb{1}X,Z\mathbb{1}Z\\},\\{XXZ,YYZ,ZZ\mathbb{1}\\},\\{XYZ,ZYX,Y\mathbb{1}Y\\}.$
(10)

Since they commute, they can be simultaneously diagonalized. The diagonalizing
unitary circuit for each cluster $V$, satisfies
$V^{\dagger}H_{p}V=\widetilde{H}_{p}$ where $p=1\cdots 5$ and
$\widetilde{H}_{p}$ is diagonal and can be written completely in terms of $Z$
and $\mathbb{1}$. Therefore, rather than doing $\exp(-iH_{p}\delta t)$, one
performs $\exp(-i(V\widetilde{H}_{p}V^{\dagger})\delta t)$. By making use of
the well-known identity:

$Z\otimes Z=\text{CNOT}\cdot(\mathbb{1}\otimes Z)\cdot\text{CNOT},$ (11)

we can change two $Z$’s to a single $Z$ and by repeatedly applying this, we
can reduce any Pauli string $\widetilde{H}_{p}$ to single $Z$. Once we have
single $Z$, then we can use the identity that $e^{itI\otimes\mathbb{P}\otimes
I}=I\otimes e^{it\mathbb{P}}\otimes I$, where $\mathbb{P}$ is a single Pauli
matrix (X, Y, or Z) because $e^{itI\otimes\mathbb{P}\otimes I}=\cos(t)I\otimes
I\otimes I+i\sin(t)I\otimes\mathbb{P}\otimes
I=I\otimes\big{(}\cos(t)I+i\sin(t)\mathbb{P}\big{)}\otimes I=I\otimes
e^{i\mathbb{P}t}\otimes I$. Therefore, $\exp(-itI\cdots IZ\cdots I\cdots
I)=I\cdots IR_{z}(2t)\cdots I\cdots I$. As an example, for the first cluster
i.e.,

$H_{1}=\alpha_{1}IZZ+\alpha_{11}ZXX+\alpha_{14}ZYY,$ (12)

we have the diagonalizing circuit made up of one CX and one Hadamard gate. It
is easy to show that:

$\widetilde{H}_{1}=\alpha_{1}IZI+\alpha_{11}ZIZ-\alpha_{14}ZZZ.$ (13)

For the first term, we do not need any CX gate, just the rotation $R_{z}$
gate. For the second term, we need two CX for the first and third qubit. For
the last term, one would naively expect four CX gates from the basic rule that
one needs $2(\nu-1)$ CX for the Pauli term with Hamming weight $\nu$, but
since the terms commute, we can use two CX gates for the second term to bring
down the required total cost to 4 CX gates. We need two CX gates for $V$ and
$V^{\dagger}$ adding to a total of 6 CX gates. The circuit for one Trotter
step is shown in Fig. 5.

Figure 5: The circuit to implement $\exp(-iH_{1}\delta t)$ for $H_{1}$ given
in (12) for one Trotter step. The CX and H gates on the left and right end are
part of the diagonalizing unitary circuit. The total cost is 6 CX gates as
mentioned in the text.

The next four clusters require 6, 4, 6, and 8 gates respectively, bringing the
total cost for $N=6$ to 30 CX gates as mentioned in Table 1 for a single
Trotter step. If we take the first cluster and do not use the information that
can be simultaneously diagonalized, it will increase the circuit complexity.
For example, if we just use the optimization levels available in QISKIT, the
cost of time evolution for (12) will be 17 CX gates while is about three times
more than our estimate. We do a similar procedure for $N=8$ and find that 70
Pauli terms can be organized into six commuting clusters each with 14, 12, 8,
14, 12, and 10 terms as given below:

$\displaystyle\\{\mathbb{1}\mathbb{1}ZZ,ZZ\mathbb{1}\mathbb{1},\mathbb{1}ZXX,Z\mathbb{1}YY,\mathbb{1}ZYY,3\mathbb{1}XX,XX\mathbb{1}Z,YYZ\mathbb{1},XXZ\mathbb{1},YY\mathbb{1}Z,XYYX,YXXY,XYXY,YXYX\\}$
$\displaystyle\\{\mathbb{1}XYZ,ZYX\mathbb{1},\mathbb{1}Y\mathbb{1}Y,ZXZX,\mathbb{1}ZYX,Z\mathbb{1}XY,X\mathbb{1}ZY,YZ\mathbb{1}X,XZXZ,Y\mathbb{1}Y\mathbb{1},XYZ\mathbb{1},YX\mathbb{1}Z\\}$
$\displaystyle\\{XYYY,YXXX,XYXX,YXYY,XXXY,YYYX,XXYX,YYXY\\}$
$\displaystyle\\{\mathbb{1}X\mathbb{1}X,ZYZY,\mathbb{1}YXZ,ZXY\mathbb{1},\mathbb{1}ZXY,Z\mathbb{1}YX,X\mathbb{1}X\mathbb{1},XZ\mathbb{1}Y,Y\mathbb{1}ZX,YZYZ,XY\mathbb{1}Z,YXZ\mathbb{1},XXXX,YYYY\\}$
$\displaystyle\\{\mathbb{1}X\mathbb{1}Y,ZYZX,\mathbb{1}Y\mathbb{1}X,ZXZY,\mathbb{1}Z\mathbb{1}Z,Z\mathbb{1}Z\mathbb{1},X\mathbb{1}Y\mathbb{1},XZYZ,Y\mathbb{1}X\mathbb{1},YZXZ,XXYY,YYXX\\}$
$\displaystyle\\{\mathbb{1}XXZ,ZYY\mathbb{1},\mathbb{1}YYZ,ZXX\mathbb{1},\mathbb{1}ZZ\mathbb{1},Z\mathbb{1}\mathbb{1}Z,X\mathbb{1}ZX,XZ\mathbb{1}X,Y\mathbb{1}ZY,YZ\mathbb{1}Y\\}.$

The quantum circuit for the time evolution of $N=8$ SYK Hamiltonian needs 22,
20, 14, 18, 16, and 20 CX gates per cluster of (Optimizing the time evolution
circuit and return probability) summing to a total of 110 CX gates per Trotter
step. After optimizing the time evolution circuits with $\delta t$, we now
need to make a good choice of the time step. For this, we did extensive tests
of the observables on noise model devices and on noiseless simulators and
found that $\delta t=1.5$ is a reasonably good choice. Our results for both
return probability and OTOC are in good agreement with this step size.

One of the observables we compute on quantum hardware is the return
probability (4). Since we focus on the infinite temperature limit ($\beta=0$)
in this work, we do not have access to the standard spectral form factor
(SFF). However, for the SYK model, certain features of the SFF are also seen
in the return probability as explored in Ref. Numasawa:2019gnl . We compute
the return probability for different $N$ using exact time evolution for
$N=6,8$. Similar to SFF, one expects to see various stages (referred to as
slope and ramp) before it converges to a constant value referred to as a
‘plateau’. The plateau value is the fraction of the initial state retained
after sufficiently long time evolution and this depends on the size of the
Hilbert space. The slope region is self-averaging (i.e., a single instance of
the model shows the same qualitative behavior as the disordered averaged case)
while the ramp and plateau regions for SFF and return probability are not
self-averaging. We show the exact results for 100 instances of the evolution
in Fig. 6. The slope region for $N=6$ ends around $t=10$ and this is
consistent with the disorder average over _only_ five instances from the
hardware results (Fig. 2) in the main text.

Figure 6: The exact result for return probability as a function of evolution
time $t$ with 100 instances of the SYK model. The initial state is taken to be
$\ket{0}^{\otimes n}$. The plateau value $\sim 1/L$ where
$L=\text{dim}(\mathcal{H})=2^{n}$ is the dimension of the Hilbert space.

## Error mitigation

The usefulness of the noisy intermediate-scale quantum (NISQ) devices is
limited by the noise and error that are associated with them. Some common
types of errors are — readout errors and gate errors. Readout errors occur
when we measure the state of the qubit incorrectly resulting in different
probabilities and expectation values. The gate errors are due to the
application of gates in a quantum circuit. They are usually classified into
coherent and incoherent errors. One useful way of differentiating between
coherent and incoherent error is how the infidelity increases with circuit
size or the number of times the channel is applied. If it increases linearly,
it is a type of incoherent noise. The coherent errors add up quadratically
2020NJPh…22g3066I . The coherent errors preserve state purity such that a
state that is lying on the surface of the Bloch sphere does not leak in the
radial direction and stays pure (on the surface). Incoherent errors map pure
states to mixed states such that they can be represented by states inside the
Bloch sphere. The incoherent errors are due to the interaction with the
outside environment. The incoherent (stochastic) errors are often modeled by
the depolarising noise model given by Urbanek:2021oej :

$\mathcal{E}(\rho)=(1-p)\rho+p\mathbb{1}/2^{n},$ (15)

where $\mathcal{E}$ denotes the completely positive trace-preserving map
(CPTP) noise channel and $\rho$ is the density matrix and $p$ is the error
rate also known as depolarising parameter. If we consider a $n$-qubit channel,
then a probability of $p=1-(1/2^{n})$ is referred to as ‘completely
depolarizing’ or ‘pure noise’ since it completely destroys the polarization of
the qubit and any pure state is changed to a maximally mixed state. The
expectation value of an operator $\mathcal{O}$ for the depolarizing channel is
given by:

$\overline{\langle O\rangle}=\text{Tr}(\mathcal{E}(\rho)O)=(1-p)\langle
O\rangle+\frac{p}{2^{n}}\text{Tr}(O)$ (16)

The dominant gate errors are coherent (do not spoil unitarity) in nature and
hence are not explained by the depolarizing noise channel. However, one can
use randomized compiling to transform them into incoherent errors. We do this
by inserting a pair of single-qubit gates on both sides such that the two-
qubit gate remains unchanged. These possible choices for the ECR gate are
shown in Table 2. The basis idea of Pauli twirling (also sometimes known as
‘randomized compiling’) He:2020udd ; PauliTwirling ; Urbanek:2021oej  is to do
operations on two-qubit gates by inserting one-qubit gates from Pauli groups
on either side (conjugation) such that the initial gate is left unchanged. The
only two-qubit gate we have used is the echoed cross-resonance (ECR) gate
which implements:

$\text{ECR}~{}q_{0},q_{1}=\frac{1}{\sqrt{2}}(IX-
XY)=\frac{1}{\sqrt{2}}\begin{pmatrix}0&0&1&-i\\\ 0&0&-i&1\\\ 1&i&0&0\\\
i&1&0&0\end{pmatrix}.$ (17)

Figure 7: The circuit implementing the twirling of the ECR gate by choosing
the $P$ gates from the Pauli group. This is also sometimes referred to as
‘randomized compiling’. This procedure changes the realistic noise which is
coherent in nature to stochastic noise.

As mentioned in the main text, there are 16 conjugations (up to a phase) that
leave the ECR gate invariant as shown in Fig. 7. One of them is the trivial
conjugation ($P_{1,2,3,4}=\mathbb{1}$). All conjugations ($P_{1}\otimes
P_{2}\cdot\text{ECR}\cdot P_{3}\otimes P_{4}$) are listed below:

$P_{1}$ | $\mathbb{1}$ | X | X | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$ | Y | Z | Z | Y | Z | Y | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathbb{1}}}$ | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$ | Y | Z | X | X
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$P_{2}$ | $\mathbb{1}$ | $\mathbb{1}$ | Z | Z | X | X | Y | Y | $\mathbb{1}$ | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$ | X | Y | Z | Z | Y | X
$P_{3}$ | $\mathbb{1}$ | X | X | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$ | Y | Z | Z | Y | Y | Z | X | X | Z | Y | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$ | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$
$P_{4}$ | $\mathbb{1}$ | $\mathbb{1}$ | Z | Z | X | X | Y | Y | Z | Z | Y | X | $\mathbb{1}$ | $\mathbb{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{1}}}$ | X | Y
Table 2: The 16 conjugations that leave the ECR gate invariant (see Fig. 7).
The conjugations in blue yield an overall phase of $\pi$.

## Protocol for OTOC computation

For the computation of OTOC as in (7), we have to make a choice of operators,
$W$ and $V$. However, the choice of operators in the four-point function does
not change the generic behavior of the growth of OTOC PhysRevLett.124.160603 .
We studied the exact time evolution using different choices of operators and
found that the qualitative difference is not substantial. Due to the limited
quantum resources in implementing the OTOC protocol and for benchmarking our
results, we have restricted to the simplest choice of $W$ and $V$ i.e., just
one non-trivial Pauli matrix in the three-qubit Pauli string.

We follow the Ref. Vermersch:2018sru  where they discussed a ‘global protocol’
to compute OTOC which uses randomly generated unitaries from the circular
unitary ensemble (CUE). By applying random unitaries, many random states are
created to mimic a thermalized scenario for the computation of the OTOC. The
protocol is developed using the following equation for traceless operators
chen2018measuring ,

$\operatorname{Tr}\left[W(t)V^{\dagger}W(t)V\right]=\frac{1}{N_{\mathcal{H}}(N_{\mathcal{H}}+1)}\overline{\langle
W(t)\rangle_{u,k_{0}}\left\langle V^{\dagger}W(t)V\right\rangle_{u,k_{0}}}.$
(18)

On the right-hand side, the overline denotes an ensemble average of
measurements over a set $u=\\{u_{0},u_{1},\cdots u_{N_{u}}\\}$ of random
unitary operators, $k_{0}$ denotes an arbitrary initial state and
$N_{\mathcal{H}}$ is the dimension of the Hilbert space. The implementation of
the global protocol requires creating a $n$-qubit random unitary operator that
applies to the input state $\ket{k_{0}}$. The decomposition of an $n$-qubit
random unitary in terms of two-qubit gates scales exponentially with $n$
Shende2006 . For a small number of qubits, the gate cost is manageable with
the current NISQ-era machine. We outline the steps to compute OTOC by the
randomized protocol developed in Ref. Vermersch:2018sru  below:

  * •

We prepare an arbitrary initial state $|k_{0}\rangle$ (position 1 in Fig.
8(a)). For simplicity, we pick the state $\ket{\mathbf{0}}=\ket{0}^{\otimes
n}$ similar to the calculation of the return probability. Next, we apply a
$n$-qubit random unitary sampled uniformly from the CUE which results in a
random state $|\psi_{1}\rangle=u\ket{\mathbf{0}}$ (position 2 in Fig. 8(a)).

  * •

The Trotter evolution of the random state with respect to the Hamiltonian is
computed using the Trotterized time-evolution operator
$U(N_{t})=\Big{[}\exp(-i\hat{H}\delta t)\Big{]}^{N_{t}}$. This yields
$|\psi_{2}\rangle=U(N_{t})|\psi_{1}\rangle$ at position 3 in the Fig. 8(a).

  * •

Finally, the necessary gates are applied to compute the observable $W$ in the
computational basis. For example, if $W=X_{i}$, an application of the Hadamard
gate before applying projective measurements on qubit $i$ in the computational
basis would allow us to compute the expectation value of $\langle
W(t)\rangle_{k_{0},u}$. The number of shots used is denoted as $N_{M}$.

  * •

Likewise, the expectation value of the operator $\langle
V^{\dagger}W(t)V\rangle_{k_{0},u}$ is computed with an inclusion of the $V$
operator for the same unitary $u$. The operator is inserted after creating the
random state $|\psi_{1}\rangle$ (position 1 in Fig. 8(b)).

  * •

The process is repeated with $N_{u}$ unitaries for each Trotter step. We
choose $N_{u}$ either to be 600 or 900 based on convergence. The unitary
matrices are drawn randomly from CUE using unitary_group.rvs available in
SciPy.

  * •

Finally, an ensemble average of the quantity $\overline{\langle
W(t)\rangle_{u,\mathbf{k}_{0}}\langle VW(t)V\rangle_{u,\mathbf{k}_{0}}}$ is
computed which is a measure of the OTOC.

We perform an additional normalization to compute the OTOC with the global
protocol

$O(t)=\frac{1}{\overline{\langle
W(t)\rangle_{u,\mathbf{k}_{0}}^{2}}}\overline{\langle
W(t)\rangle_{u,\mathbf{k}_{0}}\langle VW(t)V\rangle_{u,\mathbf{k}_{0}}}.$ (19)

It is evident that to compute the normalization there is no additional
computational cost.

Figure 8: OTOC is computed from the correlation of the measurement of two
different operators (a) $\langle W(t)\rangle$ and (b) $\langle
V^{\dagger}W(t)V\rangle$. The same set of unitaries $\\{u_{1},\cdots\
u_{N_{u}}\\}$ are required to find the correlation between the expectation
value of the operators.

Following the steps outlined above and using the measurements on circuits like
in Fig 8, we can look at the correlation of the computed expectation values.
This demonstrates the operator spreading due to information scrambling in the
system. The speed of the information spreading can thus be measured with the
computation of the OTOC using (19).

(a)

(b)

(c)

Figure 9: The change in the correlation of the operators over time in one
realization of the $N=6$ SYK model, with $\langle W(t)\rangle=\langle
Z_{1}(t)\rangle$ and $\langle VW(t)V\rangle=\langle
Z_{0}Z_{1}(t)Z_{0}\rangle$. $Z_{i}$ denotes a local Pauli-$Z$ operator on
qubit $i$.

## Additional details about simulations on IBM quantum hardware

Our return probability computation was run on ibm_nazca, ibm_cusco, and
ibm_kyoto. The OTOC computation was run on ibm_cusco and ibm_kyoto. They are
all 127 qubit devices with the Eagle r3 processor. The qubit connectivity of
these devices is shown in Fig. 10. All these machines use the same set of
basis gates which consist of ECR and single-qubit gates
$\\{\texttt{ID},\texttt{RZ},\texttt{SX},\texttt{X}\\}$.

We first construct the CNOT circuit for one Trotter step and then append the
same circuit $N_{t}-1$ times to give the circuit for $N_{t}$ Trotter steps. We
then transpile the circuits to the above basis gate set with the coupling map
of the device we used. Though, we only need $(30)N_{t}$ CNOTs for the entire
evolution as discussed in Table 1, the actual cost of the two-qubit ECR gate
is much higher. This is due to the limited connectivity of the quantum
hardware, as shown in Fig. 10. For $N=6$, we find that about 290 ECR gates are
needed for eight Trotter steps, while for $N=8$, we need about 340 ECR gates
for two Trotter steps. Thus it poses a challenge to investigate models with
non-local interactions with quantum processing units that have restricted
topology. The challenges of restricted qubit topology compared to all-to-all
qubit connectivity is discussed in the context of the multi-flavor Gross-Neveu
model in Asaduzzaman:2022bpi .

Figure 10: The general layout of the 127-qubit devices used for this work. We
shade a possible option of four utilized qubits (T-configuration) for $N=8$
which minimizes additional SWAP gates required due to the restricted
connectivity of physical qubits.

Other than minimizing SWAP gates required due to restricted topology, we use
the following criteria to map physical qubits to virtual qubits

  1. 1.

Must have low ECR gate error since two-qubit gate errors are the bottleneck in
NISQ-era devices.

  2. 2.

Have a reasonably high coherence time— for a larger number of Trotter steps,
we typically need $T_{1,2}\sim 200\,\mu$s.

  3. 3.

Preferably have a lower ECR gate application time (also sometimes referred to
as ‘gate operation time’). This enables us to perform Trotter evolution to
larger times.

The typical calibration data of some of the devices used in this work are
given in Table 3 at the time of use (November 2023).

Device |  $\text{EPLG}_{100}$ McKay:2023nxa  |  Median $T_{1}$ (in $\mu$s) |  Median $T_{2}$ (in $\mu$s) |  Median ECR gate error
---|---|---|---|---
nazca | $3.2\times 10^{-2}$ | 186.8 | 115.3 | $1.1\times 10^{-2}$
cusco | $5.9\times 10^{-2}$ | 146.6 | 83.3 | $1.74\times 10^{-2}$
kyoto | $3.6\times 10^{-2}$ | 218 | 115.4 | $8.52\times 10^{-3}$
Table 3: Comparison of some important parameters of the 127-qubit devices we
have used in this work.

Once good qubits have been chosen, we build the evolution circuits from the
Trotter circuits for a given step size. In this work, for the return
probability and the OTOC computation, we choose a Trotter step of $dt=1.5$. We
are able to use sufficiently large Trotter step size primarily because of two
reasons:

  1. 1.

The relevant quantity to consider for Trotter error estimation is
$\left(\overline{J^{2}_{ijkl}}\right)^{1/2}\,dt\sim\sqrt{3!J^{2}/N^{3}}\,dt$.
Due to the small prefactor for $J=1$, one can choose larger $dt$.

  2. 2.

The graph-coloring method reduces the number of terms in the product for
Trotterization by grouping the Pauli strings into commuting clusters, thereby
reducing the overall error, enabling larger $dt$.

For the return-probability computation, we also build a self-mitigation
circuit, where the system is evolved forward by $N_{t}/2$ steps, then evolved
back by the same number of steps. This circuit would have the same structure
as the physics circuit. We then Pauli-twirl each ECR gate in the above
circuits and generate 75 versions of the physics and self-mitigation circuits.
This makes the coherent errors in the circuit appear stochastic, and well-
described by the depolarizing channel upon averaging over all the Twirled
circuits. The mitigation circuits can be used to estimate the level of
depolarizing noise in the circuit. We run these 150 circuits for 2048 shots
each. Then, we compute the noise-less value of the observable as discussed in
(5) of the main text.