State Key Laboratory of Surface Physics and Department of Physics, Fudan
University, Shanghai 200433, China Shenzhen Institute for Quantum Science and
Engineering, Southern University of Science and Technology, Shenzhen, 518055,
China International Quantum Academy, Shenzhen, 518048, China Guangdong
Provincial Key Laboratory of Quantum Science and Engineering, Southern
University of Science and Technology, Shenzhen, 518055, China Shenzhen Key
Laboratory of Quantum Science and Engineering, Southern University of Science
and Technology, Shenzhen, 518055, China

# Quantum molecular docking with quantum-inspired algorithm

Yunting Li  Central Research Institute, Huawei Technologies, Shenzhen, 518129,
China  Xiaopeng Cui  Central Research Institute, Huawei Technologies,
Shenzhen, 518129, China  Zhaoping Xiong  Laboratory of Health Intelligence,
Huawei Cloud Computing Technologies Co., Ltd, Guizhou, 550025, China  Bowen
Liu  Central Research Institute, Huawei Technologies, Shenzhen, 518129, China
Bi-Ying Wang  Central Research Institute, Huawei Technologies, Shenzhen,
518129, China  Runqiu Shu  Central Research Institute, Huawei Technologies,
Shenzhen, 518129, China  Nan Qiao<EMAIL_ADDRESS>Laboratory of Health
Intelligence, Huawei Cloud Computing Technologies Co., Ltd, Guizhou, 550025,
China  Man-Hong Yung<EMAIL_ADDRESS>Central Research Institute,
Huawei Technologies, Shenzhen, 518129, China

###### Abstract

Molecular docking (MD) is a crucial task in drug design, which predicts the
position, orientation, and conformation of the ligand when bound to a target
protein. It can be interpreted as a combinatorial optimization problem, where
quantum annealing (QA) has shown promising advantage for solving combinatorial
optimization. In this work, we propose a novel quantum molecular docking (QMD)
approach based on QA-inspired algorithm. We construct two binary encoding
methods to efficiently discretize the degrees of freedom with exponentially
reduced number of bits and propose a smoothing filter to rescale the rugged
objective function. We propose a new quantum-inspired algorithm, hopscotch
simulated bifurcation (hSB), showing great advantage in optimizing over
extremely rugged energy landscapes. This hSB can be applied to any formulation
of objective function under binary variables. An adaptive local continuous
search is also introduced for further optimization of the discretized solution
from hSB. Concerning the stability of docking, we propose a perturbation
detection method to help ranking the candidate poses. We demonstrate our
approach on a typical dataset. QMD has shown advantages over the search-based
Autodock Vina and the deep-learning DIFFDOCK in both re-docking and self-
docking scenarios. These results indicate that quantum-inspired algorithms can
be applied to solve practical problems in the drug discovery even before
quantum hardware become mature.

##  1 Introduction

Efficient drug discovery is of significant interest in pharmaceutical science.
In the last few decades, structure-based drug design with computer-assisted
methods 1, 2, 3 has accelerated the discovery of drug candidates. _Molecular
docking_ (MD) is a crucial task in drug design, which predicts the position,
orientation, and conformation of the ligand when bound to a target protein 4,
5, 6. There are two common scenarios related to the configuration of the
ligand in MD 7: _re-docking_ , where the ligand is regarded as rigid body with
given and fixed conformation; _self-docking_ , where all degrees of freedom of
the ligand (position, orientation, and conformation) are variables. In
addition, according to the range of docking, there are also two scenarios 7:
_pocket docking_ , which is practical and widely-concerned in real experiments
with a given possible docking pocket of the protein; _blind docking_ , where
no pocket information is given and the ligand is allowed docking anywhere on
the protein.

Traditional approaches are based on systematic or stochastic search algorithms
to sample in the search space and rank the predictions according to a score
function 8. In general, the score function is usually related to potential
energy and can be physics-based, empirical, knowledge-based or machine
learning-based 9, 10. Many search-based softwares 11, 12 are available:
AutoDock 13, 14, 15, GLIDE 16, SMINA 17, and GNINA 18 et. al.. However, the
rugged energy landscape of MD problem is difficult for classical algorithms to
find the global minimum. Another method based on subgraph matching has emerged
in recent years 19, which represents both the ligand and the protein pocket as
a subgraph with pharmacophores as vertexes. According to the interaction
energy between the ligand and the protein pharmacophores, a binding
interaction graph is constructed and the MD problem is reduced to finding the
maximum weighted clique in the graph. A recent breakthrough, DIFFDOCK 20, is
based on a diffusion generative model, which defines a diffusion process
transforming the data distribution in a tractable prior and learn the score
function of the evolving distribution by machine learning. This method based
on blind docking has an impressive 30% success rate, i.e., percentage of root-
mean-square distance (RMSD) below 2Å.

In recent years, quantum computing has emerged as a promising tool to
accelerate the computation and reduce resources by exploiting quantum
mechanical systems 21, 22, 23, 24, 25, 26, 27, 28. Many attempts of structure-
based drug design with quantum computing have appeared 29, 30, 31, 32, 33, 34,
35. In the conformational generation problem, Mato et al. have applied
_quantum annealing_ (QA) 33. Several studies have focused on the protein
folding problem based on the Quantum Approximate Optimization Algorithm (QAOA)
36, 37, 38, 39, 40. In addition, Pandey et al. have realized a simplized
multibody docking problem on a quantum annealer 34. To address molecular
docking problem, we explore the effectiveness of quantum-inspired algorithms
41, 23, 24, 42 for support, which are expected to tunnel through ”high but
narrow” barriers in energy landscapes. Specifically, we employ the simulated
bifurcation (SB) algorithm, which is a purely quantum-inspired algorithm based
on quantum adiabatic optimization using nonlinear oscillators 43, 44, 45, 46,
47, 48. SB is driven by Hamiltonian dynamics evolution, and allows for the
simultaneous updating of variables, resulting in great parallelizability and
accelerated combinatorial optimization. SB also shows significant potential
for finding the global minimum. Recently, two SB variants, namely _ballistic
SB_ (bSB) and _discrete SB_ (dSB), have been proposed for further enhancement
46.

In this work, we propose an efficient _quantum molecular docking_ (QMD)
approach to solve MD problem with quantum-inspired algorithms. The Gray code
method is applied to discrete and encode the degrees of freedom of ligand into
binary variables. We also construct a phase encoding method to efficiently
encode angular variables. We determine the interaction energy function as the
objective function in MD optimization problem, and propose a smoothing filter
to rescale the rugged objective function, which significantly benefits further
optimization. In addition, we propose a new SB variant, _hopscotch SB_ (hSB),
which has great advantage in optimizing over rugged energy landscapes. This
hSB can be applied to any formulation of objective function under binary
variables. An adaptive local continuous search is also introduced for further
optimization of the discretized solution from hSB. Concerning the stability of
docking, we propose a perturbation detection method to help ranking the
candidate poses. We demonstrate the performance of our QMD approach on the
same testing dataset as DIFFDOCK with respect to RMSD and success rate
metrics. In pocket re-docking scenario, QMD increases the top-1 success rate
of the root-mean-square distance (RMSD) less than 2Å by about 6% compared with
DIFFDOCK. In pocket self-docking scenario, QMD with perturbation detection
increases the top-1 success rate by about 2% compared with Autodock Vina, and
increases the top-5 success rate of RMSD less than 5Å by about 6% compared
with DIFFDOCK.

##  2 Background

###  2.1 Molecular docking problem

Molecular docking problem aims to predict the position, orientation and
conformation of a ligand (usually referred to as pose) when bound to a target
protein. A score function is introduced to estimate the correctness of an
optimized pose and rank several optimized samples.

In general, the target protein is rigid and the structure is given. The self-
docking allows the degrees of freedom of the ligand including three
translations $T_{x},T_{y},T_{z}$ in Cartesian coordinate system, three
rotation angles $\phi_{x},\phi_{y},\phi_{z}$ around the ligand itself, and $M$
torsion angles $\theta_{1},\dots,\theta_{M}$ along $M$ rotatable bonds of the
ligand. The re-docking refers to only consider the first two rigid
transformations (translation and rotation) with a fixed conformation given as
input.

Suppose the position of an atom $i$ in the ligand is denoted by $\vec{r}_{i}$.
Without loss of generality, the ligand is updated in the sequence of
rotations, translations and torsions. The translations and rotations update
the atom position as follows:

$\displaystyle\vec{r}^{\prime}_{i}=R(\phi_{x},\phi_{y},\phi_{z})(\vec{r}_{i}-\vec{r}_{0})+\vec{r}_{0}+\vec{T},$
(1)

where $\vec{r}_{0}$ is the geometric center of the current ligand,
$\vec{T}=(T_{x},T_{y},T_{z})^{T}$ is a vector giving the translation of the
origin, and
$R(\phi_{x},\phi_{y},\phi_{z})=R_{z}(\phi_{z})R_{y}(\phi_{y})R_{x}(\phi_{x})$
is the rotation matrix obtained from three basic rotation matrices of angles
$\phi_{x},\phi_{y},\phi_{z}$ around $x,y,z$ axis. Next the torsions update the
atom position $\vec{r}^{\prime}_{i}$ as follows:

$\displaystyle\vec{r}^{\prime\prime}_{i}=\mathbf{R}_{M}\circ\cdots\circ\mathbf{R}_{2}\circ\mathbf{R}_{1}(\vec{r}^{\prime}_{i}),$
(2)

where $\mathbf{R}_{i},i=1,\dots,M$ is the torsion update of angle $\theta_{i}$
along the $i$th rotatable bond. The details of torsion update are elaborated
in Sec. S1 in the Supplemental Material 49.

The molecular docking process can be viewed as an optimization problem:
minimizing the interaction energy over the degrees of freedom of the ligand to
generate a stable complex of ligand and protein. The interaction energy
between the ligand and protein includes the Lennard-Jones potential energy and
the electrostatic potential energy, which is usually regarded as the objective
function and the score function.

To be specific, the Lennard-Jones potential models soft repulsive and
attractive (van der Waals) interactions for electronically neutral atoms or
molecules. The commonly used expression reads 14, 15

$\displaystyle
V_{\mathrm{LJ}}=\sum_{i,j}V_{\mathrm{LJ}}^{ij}=\sum_{i,j}\varepsilon_{ij}(\frac{R_{ij}}{r_{ij}})^{12}-2\varepsilon_{ij}(\frac{R_{ij}}{r_{ij}})^{6},$
(3)

where $i,j$ represent the atoms in the ligand and the protein respectively,
$R_{ij}$ is the sum of radius of atom $i$ and $j$, $r_{ij}$ is the distance
between atom $i$ and $j$, and $\varepsilon_{ij}$ is the depth of the potential
well of atom $i$ and $j$.

The electrostatic potential energy results from conservative Coulomb forces
and is associated with the configuration of a particular set of point charges
within a defined system. The commonly used expression reads 14, 15

$\displaystyle
V_{e}=\sum_{i,j}V_{e}^{ij}=\sum_{i,j}\frac{q_{i}q_{j}}{f(r_{ij})r_{ij}},$ (4)

where $q_{i},q_{j}$ are the charges of atom $i$ and $j$ respectively, $r_{ij}$
is the distance between atom $i$ and $j$, and $f(r_{ij})$ is a correction
function of $r_{ij}$.

Therefore, the total interaction energy is a weighted intermolecular pair
potential reads 14, 15

$\displaystyle V=w_{1}V_{\mathrm{LJ}}+w_{2}V_{e},$ (5)

where $w_{1},w_{2}$ are the weights for two potential energy. The degrees of
freedom of the ligand can affect positions of each atom resulting in the
change of $r_{ij}$ terms. The optimization problem is to minimize the total
interaction energy $V$ as the objective function. The details of parameters,
the correction function and other techniques are elaborated in Sec. S2 in the
Supplemental Material 49.

###  2.2 Quantum-inspired algorithms

The original simulated bifurcation (SB) 43, 44, 45, 46 is a heuristic
technique for accelerating combinatorial optimization in the Ising
formulation. A network of quantum-mechanical nonlinear oscillators _Kerr-
nonlinear parametric oscillators_ (KPOs) is introduced in the quantum
Hamiltonian 43, 44. Driven by the derived classical Hamiltonian dynamics
evolution from Hamilton-Jacobi equation, the equations of motion of position
$x_{i}$ and momentum $y_{i}$ of the $i$th particle can be obtained. All
positions and momentums are randomly set around zero at the initial time.
During the evolution, the positions $x_{i}$ of each particle oscillate and
then move towards values $\pm 1$, exhibiting as bifurcations. Finally, the
$\operatorname{sgn}(x_{i})$ gives the solution of the corresponding spin value
$s_{i}$ in the Ising formulation.

From the intuition of binary spin values, a variant bSB is introduced 46 with
inelastic walls set at $x_{i}=\pm 1$. Moreover, the discrete version of bSB,
namely, dSB is introduced 46 to suppress analog deviation by discretizing
$x_{i}$ with $\operatorname{sgn}(x_{i})$ when updating the derivative of the
momentum $y_{i}$.

Fortunately, recent study has shown that dSB can be potentially generalized to
higher-order unconstrained binary objective (HUBO) formulation and
demonstrated the advantage in a three-order task 47. By replacing the Ising
energy term with a general higher-order objective function $E(\vec{x})$, we
have the general equations of motion as follows:

$\displaystyle\dot{x}_{i}$ $\displaystyle=\frac{\partial
H_{\mathrm{dSB}}}{\partial y_{i}}=a_{0}y_{i}$ (6) $\displaystyle\dot{y}_{i}$
$\displaystyle=-\frac{\partial H_{\mathrm{dSB}}}{\partial
x_{i}}=\left[a_{0}-a(t)\right]x_{i}+c_{0}\left.\frac{\partial
E(\vec{x})}{\partial x_{i}}\right|_{\vec{x}=\operatorname{sgn}(\vec{x})},$

where $\vec{x}$ is the vector of positions of all spins, $x_{i},y_{i}$ are the
positions and momentums of a particle corresponding to the $i$th spin, $a(t)$
is a control parameter increased from zero to $a_{0}$, and $a_{0},c_{0}$ are
positive constants. The Hamiltonian $H_{\mathrm{dSB}}$ reads

$H_{\mathrm{dSB}}=\frac{a_{0}}{2}\sum_{i}^{N}y_{i}^{2}+V_{\mathrm{dSB}},$ (7)

and if $|x_{i}|\leq 1$ for all $x_{i}$, the potential $V_{\mathrm{dSB}}$ reads

$V_{\mathrm{dSB}}=\frac{a_{0}-a(t)}{2}\sum_{i}^{N}x_{i}^{2}-c_{0}E(\vec{x}),$
(8)

otherwise it reads infinity.

The algorithmic advantages of SB are that the equation of motion is easy to
implement with only product-sum operation given the derivatives; the variables
can be updated simultaneously for parallel acceleration; and SB does not
contain probabilistic processes, and therefore can be easily hardwired 44.

Since SB algorithms require binary variables, encoding methods, including
binary encoding, Gray encoding 50 and one-hot encoding 51, are necessary to
encode general variables into binary variables in general tasks. A suitable
choice of encoding methods may also improve the efficiency and accuracy for
the particular task.

##  3 Methods

###  3.1 Overview

Figure 1: The framework of Quantum-inspired Molecular Docking. The inputs are
the 3D representations of the ligand and protein which can be generated by an
optional module (molecular conformational generation). The degrees of freedom
are discretized and encoded to the binary variable formulation. The objective
function is of binary variables based on potential energy and then transformed
through a smoothing filter. After optimized by quantum-inspired algorithms,
the optimized binary result is decoded to discrete values of the degrees of
freedom. The local continuous search is an optional module for further
optimization in the continuous domain. The optimized continuous result is
output to recover the optimized pose of the ligand. The perturbation detection
is an optional module for helping rank several optimized poses by virtue of
the stability.

The framework in Figure 1 shows a general process of quantum-inspired
molecular docking. The inputs are the 3D representations of the ligand and the
protein. If only an 1D/2D representation or a rough 3D representation is
available, we introduce an optional module, molecular conformational
generation, to generate the 3D representation based on potential energy 52 or
molecular volume 33. The degrees of freedom of the ligand including
translations $T_{x},T_{y},T_{z}$, rotations $\phi_{x},\phi_{y},\phi_{z}$ and
torsions $\theta_{1},\dots,\theta_{M}$ are discretized and encoded to binary
variables. The objective function is of binary variables based on potential
energy and then transformed through a smoothing filter. After optimized by
quantum-inspired algorithms in the binary domain, the solution is then decoded
to the corresponding discrete values of the degrees of freedom. For further
optimization in the continuous domain, an optimized continuous solution is
given by a local continuous search. Finally, the output recovers the optimized
pose of the ligand.

After we obtain several optimized poses, the best pose is ranked by the score
function. Considering the stability of the crystallization in docking
experiments, the perturbation detection is an optional module to better rank
the candidate poses.

###  3.2 Discretization formulation

The degrees of freedom can be discretized into uniformly discrete values by
encoding methods. Here we generalize the gray code method with coefficients to
encode any range of uniformly discrete values. We also adapt the phase
encoding method to encode trigonometric functions. If the degree of freedom
corresponds to $n$ binary variables, then it has ${d=2^{n}}$ uniformly
discrete values. The gray code is suitable for all three kinds of variables
(translation, rotation and torsion), while phase encoding is suitable for
encoding angle variables (rotation and torsion).

####  3.2.1 Gray code method

In gray code method 50, the bit strings of two adjacent integers differ in
only one bit. Given $n$ binary variables
$\mathbf{s}=(s_{0},\dots,s_{n-1})\in\\{-1,1\\}^{n}$, the correspondence to
${d=2^{n}}$ integers ${\\{-d/2+1,-d/2+2,\dots,d/2\\}}$ is as follows:

$\displaystyle
D(\mathbf{s})=\frac{1}{2}\left(\sum_{j=0}^{n-1}2^{n-j-1}\prod_{i=0}^{j}s_{i}+1\right).$
(9)

Suppose the step sizes of translation, rotation and torsion are $\Delta
T,\Delta\phi,\Delta\theta$ respectively. Then the discrete values of each
degree of freedom are $\Delta TD(\mathbf{s}),\Delta\phi
D(\mathbf{s}),\Delta\theta D(\mathbf{s})$. To ensure sufficient sampling,
$\Delta T$ is at most 0.5Å and $\Delta\phi,\Delta\theta$ are about
$5^{\circ}$. Here we set $\Delta T=T/d$ where $T$ is the translation range,
and $\Delta\phi=\Delta\theta=2\pi/d$. In general, the value of $d$ can be
different for translation, rotation and torsion depending on the needs.

####  3.2.2 Phase encoding method

In phase encoding method 53, the correspondence between the phase
$\mathrm{e}^{\mathrm{i}\psi}$ of an angle $\psi\in\\{\psi_{k}\\}$ and a
polynomial $p_{n}(\mathbf{s})$ is established. The detail is shown in the
following lemma from Ref. 53.

######  Lemma 1.

Suppose $n$ is a positive integer and $\\{\psi_{k}\\}_{k=0}^{2^{n}-1}$ are
phases uniformly distributed on the circle. Then we can construct a polynomial
$p_{n}(\mathbf{s})$ with $2^{n-1}$ terms such that
$\\{p_{n}(\mathbf{s})|\mathbf{s}\in\\{\pm
1\\}^{n}\\}=\\{\mathrm{e}^{\mathrm{i}\psi_{k}}\\}_{k=0}^{2^{n}-1}$.

Since the angles only appear in the rotation operator in the form of
$\sin\psi,\cos\psi$, then we propose the efficient polynomial expressions of
$\sin\psi,\cos\psi$ from the polynomial expression $p_{n}(\mathbf{s})$ of
$\mathrm{e}^{\mathrm{i}\psi}$ read

$\displaystyle\sin\psi=\Im p_{n}(\mathbf{s}),\quad\cos\psi=\Re
p_{n}(\mathbf{s}).$ (10)

In addition, we can also obtain the efficient polynomial expression of $\psi$
as follows:

$\displaystyle\psi=-\mathrm{i}\ln p_{n}(\mathbf{s}).$ (11)

If $n=3$ for example, we have $\psi_{k}=\frac{2\pi}{2^{n}}k=\frac{\pi}{4}k$
and the corresponding polynomial reads

$p_{3}(\mathbf{s})=c_{0}s_{0}+c_{1}s_{1}+c_{2}s_{2}+c_{3}s_{0}s_{1}s_{2},$
(12)

where $\mathbf{s}=(s_{0},s_{1},s_{2})\in\\{\pm 1\\}^{3}$ and the coefficients
are

$\begin{pmatrix}c_{0}\\\ c_{1}\\\ c_{2}\\\
c_{3}\end{pmatrix}=\frac{1}{4}\begin{pmatrix}1+(\sqrt{2}-1)\mathrm{i}\\\
1-(\sqrt{2}+1)\mathrm{i}\\\ (1+\sqrt{2})+\mathrm{i}\\\
(1-\sqrt{2})+\mathrm{i}\end{pmatrix}.$ (13)

Then the polynomial expressions of $\sin\psi,\cos\psi$ can be obtained from
straightforward calculations. Similarly, for any value of $n$, we can obtain
the polynomial expression of discrete phases according to Ref. 53.

Therefore, replacing the denotation $\psi$ with the concerning rotation angles
$\phi_{x},\phi_{y},\phi_{z}$ and torsion angles $\theta_{1},\dots,\theta_{M}$,
we obtain the encoding correspondence between the $d$ discrete values and the
$n$ binary variables $\mathbf{s}$. The step sizes of rotation angles and
torsion angles are $\Delta\phi=\Delta\theta=2\pi/d$. In general, the value of
$d$ can also be different for rotation and torsion depending on the needs.

Figure 2: The effect of the smoothing filter. We set $V_{0}=10$, $H_{0}=-10$,
$\Delta V=20$ and thus $x_{0}=10$. The two dashed vertical lines separate the
domain into three areas corresponding to the three cases in Eq. (14)
respectively.

###  3.3 Smoothing filter

To further facilitate the optimization of rough objective function, we propose
the smoothing filter technique to compress the huge peak and expand the sharp
minimum. Since in MD the objective function is inversely proportional to the
12th power of the distance between two atoms in Eq. (5), the close atom pair
will cause huge value of objective function, which is not suitable for
traditional optimization algorithms.

The smoothing filter is realized by a map function. We map the original
objective function $V$ defined in Eq. (5) to a bounded range
$(-V_{0}-1,V_{0}+1)$, $V_{0}>0$. The map function is designed to maintain the
monotonicity, flatten the objective function around a given low enough value,
and greatly compress the huge peak of the objective function. The form of map
functions is not unique and can be constructed differently according to the
specific problem.

Here we construct a map function as follows:

$\displaystyle
V^{\prime}=\begin{cases}\exp[\frac{\log(V_{0}+1)}{x_{0}-H_{0}}(V-H_{0})]-(V_{0}+1),\quad\mathrm{if}\
V\leq x_{0},\\\ \log V-\log x_{0}+x_{0},\quad\mathrm{if}\ x_{0}<V\leq 10^{\log
x_{0}+V_{0}},\\\ \tanh(\log V-\log x_{0}-V_{0})+V_{0},\quad\mathrm{if}\
V>10^{\log x_{0}+V_{0}},\end{cases}$ (14)

where $H_{0}$ is the given low enough value with the similar magnitude of the
minimum, $x_{0}={\max\\{H_{0}+\Delta V,1\\}}$ is the upper bound for
flattening, and $\Delta V$ is the range of flattening. Once $V\geq x_{0}$, $V$
usually have huge values scaling as power of 10 due to the overlap of atom
pairs. The effect of the smoothing filter is shown in Figure 2. Here we set
$V_{0}=10$, $H_{0}=-10$, $\Delta V=20$ and $x_{0}=10$ for illustration, which
is also suitable for molecular docking problem. The two dashed vertical lines
separate the domain into three areas corresponding to the three cases in Eq.
(14) respectively.

Therefore, the objective function is transformed through the smoothing filter
and then is used for optimization.

###  3.4 Hopscotch simulated bifurcation

In this work, we propose a new variant of SB, hSB, for optimizing rough
objective function with sharp minimums and huge peaks. Additionally, hSB does
not require the objective function to be in HUBO formulation and does not
require the information of derivatives, and therefore is capable for any
formulation under binary variables.

Even though the partial derivative of $E(\vec{x})$ is calculated at
$\vec{x}=\operatorname{sgn}(\vec{x})$ in dSB as shown in Eq. (6), it still
comes from the continuous domain. This may give wrong evolution directions
since the derivative of a continuous position is not in correspondence with
the one of discrete difference as shown in Figure 3. The dark red dot (the
most left) represents the global minimum. When evolving to the red circle
position, the derivative of the continuous objective function will indicate
the black arrow as the next evolution direction. However, the discrete
difference calculated from the nearest two red dots will indicate the red
arrow which is correct. When there is a peak between two discrete positions,
bSB and dSB may lead to wrong evolution path and miss the global minimum.

To avoid the misleading of information in the continuous domain, the evolution
of hSB only depends on the discrete positions $\operatorname{sgn}(x_{i})$. We
replace the partial derivative $\left.\frac{\partial E(\vec{x})}{\partial
x_{i}}\right|_{\vec{x}=\operatorname{sgn}(\vec{x})}$ with the difference of
adjacent discrete positions as follows:

$\displaystyle-\frac{1}{2}\left(E[e_{i}(\operatorname{sgn}(\vec{x}))]-E[\operatorname{sgn}(\vec{x})]\right)\operatorname{sgn}(x_{i}),$
(15)

where $e_{i}(\operatorname{sgn}(\vec{x}))$ is $\operatorname{sgn}(\vec{x})$
with the $i$th element flipped. Therefore, the influence from the non-discrete
region is completely removed, and the optimization strictly depends on the
objective function in discrete region.

The hSB yields great improvements in optimizing rough objective functions and
allows any formulation of the objective function. Additionally, all variables
can be simultaneously updated at each iteration step if parallel computation
is used. The time cost of hSB can be theoretically reduced by the times of the
number of variables if parallel computation is available.

Figure 3: The schematic diagram of the improvement in hSB. The variable space
is discretized. The dark red dot (the most left) represents the global
minimum. The black arrow indicates the derivative of the continuous objective
function at the red circle position, while the red arrow indicates the
discrete difference calculated from the nearest two red dots.

###  3.5 Local continuous search

After the discrete optimization by quantum-inspired algorithm hSB, we propose
a local continuous search to further optimize over the degrees of freedom of
translations $T_{x},T_{y},T_{z}$, rotations $\phi_{x},\phi_{y},\phi_{z}$ and
torsions $\theta_{1},\dots,\theta_{M}$ in the continuous domain.

Here we consider the adaptive gradient descent (Adagrad) with adaptive
learning rate. Denote by $f(\Theta)$ the objective function and $\Theta$ the
variables. The evolution of Adagrad reads

$\displaystyle G_{t}$ $\displaystyle=G_{t-1}+\nabla f(\Theta),$ (16)
$\displaystyle\Theta_{t+1}$
$\displaystyle=\Theta_{t}-\frac{\alpha_{t}}{\sqrt{G_{t}+\epsilon}}\nabla
f(\Theta),$

where $G_{t}$ is a vector of the squared sum of the past-time gradients with
regards to each variable, $\epsilon$ is a correction term to avoid dividing by
0 and is generally insignificantly small ($10^{-8}$), and $\alpha_{t}$ is an
adaptive learning rate to adjust the weight of $G_{t}$ for accelerating the
optimization.

###  3.6 Perturbation detection

In general, the interaction energy is usually used as the score function to
rank several optimized poses. However, the crystal pose is not exactly
corresponding to the lowest-energy pose. In the docking experiments, the
crystallization process is also related to the stability of the pose.
Therefore, we propose the perturbation detection to help rank the optimized
poses.

For each optimized pose, we slightly perturb the pose for $M$ times by adding
a small translation and rotation and record the corresponding interaction
energy after each perturbation. We calculate the standard deviation of the $M$
interaction energy. When ranking the optimized poses, we first filter the $k$
lowest energy within a small $\delta$ difference (e.g. $\delta=5\%$), then the
best pose is the one with the smallest standard deviation within the $k$
poses. Note that if the lowest energy is greatly smaller than others, i.e.,
the difference is larger than $\delta$, then the best pose is still the one
with the lowest energy.

##  4 Results

###  4.1 Re-docking scenario

In re-docking scenario, the target protein and the ligand are regarded as
rigid bodies. Thus the ligand has three degrees of freedom of translation and
three degrees of freedom of rotation. The ligand conformation is pre-generated
and given as the input. Here we compare our approach QMD with the ML-based
method DIFFDOCK. During the inference process in DIFFDOCK, the ligand
configuration can be fixed. We evaluate on the same testing dataset of 363
protein-ligand pairs as DIFFDOCK 20. For re-docking, the approaches receive
the crystal structures of the protein and the ligand as inputs. 10 initial
poses of the ligand are randomly selected, and the corresponding predicted
docking pose are ranked according to the score function in Eq. (5). The
implementation details of QMD and Autodock Vina are explained in Sec. S3 in
the Supplemental Material 49.

Table 1: Re-docking comparison between QMD and DIFFDOCK. The top-1 RMSD is the RMSD of the highest ranked pose; the top-5 RMSD is the smallest RMSD among the 5 highest ranked poses. Shown is the percentage of predicted poses with RMSD $<2$Å and the mean RMSD for top-1 and top-5 respectively. Method | Top-1 RMSD (Å) | Top-5 RMSD (Å)
---|---|---
%$<2$ | Mean | %$<2$ | Mean
DIFFDOCK (10) | $49.94\pm 0.78$ | $7.66\pm 0.21$ | $59.50\pm 1.40$ | $3.99\pm 0.05$
QMD (10, no local) | $46.40\pm 2.65$ | $4.16\pm 0.18$ | $56.11\pm 2.63$ | $2.64\pm 0.11$
QMD (10) | $55.37\pm 1.68$ | $3.59\pm 0.04$ | $63.77\pm 2.04$ | $2.37\pm 0.09$

We consider the RMSD between the predicted docking pose and the ground-truth
(crystal pose). The comparison of re-docking is shown in Table 1. The top-1
RMSD is the RMSD of the highest ranked pose; the top-5 RMSD is the smallest
RMSD among the 5 highest ranked poses, which is a useful metric when multiple
predictions are used for downstream verification and applications. The success
rate (i.e., percentage of predicted poses with RMSD $<2$Å) and the mean RMSD
for top-1 and top-5 are shown respectively. The experiments are repeated for 3
times and the average success rate with standard deviation is shown. QMD
improves the success rate about 5% compared to DIFFDOCK, and the local
continuous search in QMD can improve the success rate about 9% based on the
results from hSB.

###  4.2 Self-docking scenario

We conduct a comparative analysis between our approach with the search-based
method Autodock Vina and DIFFDOCK in the pocket self-docking scenario, where
the ligand has the degrees of freedom of translation, rotation and torsion.
Autodock Vina is implemented by the open Python package, ODDT (Open Drug
Discovery Toolkit) 54. We evaluate on the same testing dataset of 363 protein-
ligand pairs as DIFFDOCK 20. For self-docking, the approaches receive the
crystal pose of the protein and the generated conformation of the ligand by
UFF (Universal Force Field) optimization in RDKit. The pocket is of size
$12\times 12\times 12$Å3 centered at the crystal pose. For each approach, 10
initial poses of the ligand are randomly selected, and the corresponding
predicted docking pose are ranked according to the score function. In
addition, we also consider the QMD with perturbation detection. The
implementation details are explained in Sec. S3 in the Supplemental Material
49.

Table 2: Self-docking comparison between QMD and Autodock Vina. The top-1 RMSD is the RMSD of the highest ranked pose; the top-5 RMSD is the smallest RMSD among the 5 highest ranked poses. Shown is the percentage of predicted poses with RMSD $<2$Å and the mean RMSD for top-1 and top-5 respectively.  Method | Top-1 RMSD (Å) | Top-5 RMSD (Å)
---|---|---
% $<2$ | Mean | % $<2$ | Mean
DIFFDOCK (10) | $30.53\pm 0.79$ | $7.58\pm 0.32$ | $38.56\pm 1.10$ | $4.81\pm 0.21$
Autodock Vina (10) | $27.36\pm 0.94$ | $4.83\pm 0.04^{*}$ | $46.12\pm 1.71$ | $3.40\pm 0.20^{*}$
QMD (10, no local) | $23.58\pm 0.84$ | $5.74\pm 0.11$ | $37.88\pm 0.71$ | $3.59\pm 0.04$
QMD (10) | $27.93\pm 1.14$ | $5.46\pm 0.04$ | $43.18\pm 1.45$ | $3.28\pm 0.05$
QMD (10, perturbation) | $29.36\pm 0.81$ | $5.38\pm 0.10$ | $45.45\pm 1.98$ | $2.92\pm 0.02$

We consider the RMSD between the predicted docking pose and the ground-truth
(crystal pose). The comparison of self-docking is shown in Table 2. The top-1
RMSD is the RMSD of the highest ranked pose; the top-5 RMSD is the smallest
RMSD among the 5 highest ranked poses. We show the success rate (i.e.,
percentage of predicted poses with RMSD $<2$Å) and the mean RMSD for top-1 and
top-5 respectively. The experiments are repeated for 3 times and the average
success rate with standard deviation is shown. QMD with perturbation detection
improves the top-1 success rate about 2% compared to Autodock Vina, and
improves the top-5 success rate about 6% compared to DIFFDOCK. The local
continuous search in QMD can improve the top-1 success rate about 4% based on
the results from hSB, and the perturbation detection can improve the success
rate about 2%.

However, the small standard deviation of the mean RMSD of Autodock Vina is
because for several samples, Autodock Vina can not output any valid docking
result with its embedding selection process. Therefore, the mean RMSD is
averaged over obtained results usually with good RMSDs, which is not the mean
RMSD for the whole dataset.

The fact that the crystal pose is not guaranteed to be the lowest-energy pose
under such a model explains the difficulty to optimize the docking pose. We
find that among the 363 testing protein-ligand pairs, about 65% predicted
docking poses in QMD have lower interaction energy than the crystal poses.
This implies that our approach has the ability to find the solution with as
low as possible objective function, which unfortunately is not the crystal
pose. The fine construction of the score function may be an independent and
important research in the structure-based drug design, but is not the main
focus in our study.

Further, compared with the re-docking scenario shown in Table 1, the higher
success rate is due to the rigid-body condition which constrains the number of
possible poses. However, in self-docking with the degree of freedom of
torsion, it may find other better solutions by virtue of a more suitable
configuration, which can not be reached in re-docking.

##  5 Conclusions

In this work, we proposed an efficient QMD approach to solve MD problem. We
constructed two binary encoding methods to efficiently discretize the degrees
of freedom. We determined the interaction energy function as the objective
function, and proposed a smoothing filter to rescale the rugged objective
function, which significantly benefits further optimization. We proposed a new
hSB algorithm showing great advantage in optimizing over rugged energy
landscapes. An adaptive local continuous search is introduced for further
optimization of the discretized solution from hSB. We also proposed the
perturbation detection to help rank the optimized poses, which considers the
stability of docking and effectively improves the success rate. We
demonstrated the performance of our QMD approach on the same testing dataset
as DIFFDOCK with respect to RMSD and success rate metrics. In pocket re-
docking scenario, QMD increases the top-1 success rate of the root-mean-square
distance (RMSD) less than 2Å by about 6% compared with DIFFDOCK. In pocket
self-docking scenario, QMD with perturbation detection increases the top-1
success rate by about 2% compared with Autodock Vina, and increases the top-5
success rate of RMSD less than 5Å by about 6% compared with DIFFDOCK.

This work indicates that quantum-inspired algorithms have potentials to solve
practical problems even before quantum hardware become mature, including
molecular configuration conformation, molecular docking and protein folding
problems in the drug discovery field. More practical scenarios have appeared
to seek solution possibilities of QAOA, quantum annealing, and quantum
computation in the near future. This work opens a new direction for solving MD
problem. In the future, it is also worth studying the generalization of QMD in
blind docking scenario. Customized parameter settings based on the
characteristics of the target ligand and protein can also be considered.

Many technique details including the binary encoding methods, the smoothing
filter, and the hSB algorithm in our work have numerous applications in other
similar scenarios. The generalization of hSB in other classes of problems is
worth further study.

## Acknowledgments

Y. Li is grateful to S. Pan and J. Hu for insightful discussions. B. Liu is
partially supported by the National Natural Science Foundation of China
(Grants No. 12101394 and No. 12171426)

## References

  * Leelananda and Lindert 2016 Leelananda, S. P.; Lindert, S. Computational methods in drug discovery. _Beilstein J. Org. Chem._ 2016, _12_ , 2694–2718, PMC5238551
  * dos Santos et al. 2018 dos Santos, R. N.; Ferreira, L. G.; Andricopulo, A. D. In _Computational Drug Discovery and Design_ ; Gore, M., Jagtap, U. B., Eds.; Springer New York: New York, NY, 2018; pp 31–50
  * Mulligan 2020 Mulligan, V. K. The emerging role of computational design in peptide macrocycle drug discovery. _Expert Opin. Drug Discov._ 2020, _15_ , 833–852, PMID: 32345066
  * Kuntz et al. 1982 Kuntz, I. D.; Blaney, J. M.; Oatley, S. J.; Langridge, R.; Ferrin, T. E. A geometric approach to macromolecule-ligand interactions. _J. Mol. Biol._ 1982, _161_ , 269–288
  * Gschwend et al. 1996 Gschwend, D. A.; Good, A. C.; Kuntz, I. D. Molecular docking towards drug discovery. _J. Mol. Recognit._ 1996, _9_ , 175–186
  * Huang et al. 2006 Huang, N.; Shoichet, B. K.; Irwin, J. J. Benchmarking Sets for Molecular Docking. _J. Med. Chem._ 2006, _49_ , 6789–6801, PMID: 17154509
  * Lu et al. 2022 Lu, W.; Wu, Q.; Zhang, J.; Rao, J.; Li, C.; Zheng, S. TANKBind: Trigonometry-Aware Neural NetworKs for Drug-Protein Binding Structure Prediction. _NeurIPS_ 2022,
  * Stanzione et al. 2021 Stanzione, F.; Giangreco, I.; Cole, J. C. In _Chapter Four - Use of molecular docking computational tools in drug discovery_ ; Witty, D. R., Cox, B., Eds.; Progress in Medicinal Chemistry; Elsevier, 2021; Vol. 60; pp 273–343
  * Li et al. 2014 Li, Y.; Han, L.; Liu, Z.; Wang, R. Comparative Assessment of Scoring Functions on an Updated Benchmark: 2. Evaluation Methods and General Results. _J. Chem. Inf. Model._ 2014, _54_ , 1717–1736
  * Li et al. 2019 Li, J.; Fu, A.; Zhang, L. An Overview of Scoring Functions Used for Protein–Ligand Interactions in Molecular Docking. _Interdiscip. Sci._ 2019, _11_ , 320–328
  * Wang et al. 2016 Wang, Z.; Sun, H.; Yao, X.; Li, D.; Xu, L.; Li, Y.; Tian, S.; Hou, T. Comprehensive evaluation of ten docking programs on a diverse set of protein–ligand complexes: the prediction accuracy of sampling power and scoring power. _Phys. Chem. Chem. Phys._ 2016, _18_ , 12964–12975
  * Pagadala et al. 2017 Pagadala, N. S.; Syed, K.; Tuszynski, J. Software for molecular docking: a review. _Biophys. Rev._ 2017, _9_ , 91–102
  * Morris et al. 1998 Morris, G. M.; Goodsell, D. S.; Halliday, R. S.; Huey, R.; Hart, W. E.; Belew, R. K.; Olson, A. J. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. _J. Comput. Chem._ 1998, _19_ , 1639–1662
  * Trott and Olson 2010 Trott, O.; Olson, A. J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. _J. Comput. Chem._ 2010, _31_ , 455–461
  * Eberhardt et al. 2021 Eberhardt, J.; Santos-Martins, D.; Tillack, A. F.; Forli, S. AutoDock Vina 1.2.0: New Docking Methods, Expanded Force Field, and Python Bindings. _J. Chem. Inf. Model._ 2021, _61_ , 3891–3898
  * Friesner et al. 2004 Friesner, R. A.; Banks, J. L.; Murphy, R. B.; Halgren, T. A.; Klicic, J. J.; Mainz, D. T.; Repasky, M. P.; Knoll, E. H.; Shelley, M.; Perry, J. K.; Shaw, D. E.; Francis, P.; Shenkin, P. S. Glide: A New Approach for Rapid, Accurate Docking and Scoring. 1. Method and Assessment of Docking Accuracy. _J. Med. Chem._ 2004, _47_ , 1739–1749
  * Koes et al. 2013 Koes, D. R.; Baumgartner, M. P.; Camacho, C. J. Lessons Learned in Empirical Scoring with smina from the CSAR 2011 Benchmarking Exercise. _J. Chem. Inf. Model._ 2013, _53_ , 1893–1904
  * McNutt et al. 2021 McNutt, A. T.; Francoeur, P.; Aggarwal, R.; Masuda, T.; Meli, R.; Ragoza, M.; Sunseri, J.; Koes, D. R. GNINA 1.0: molecular docking with deep learning. _J. Cheminformatics_ 2021, _13_ , 43
  * Banchi et al. 2020 Banchi, L.; Fingerhuth, M.; Babej, T.; Ing, C.; Arrazola, J. M. Molecular docking with Gaussian Boson Sampling. _Sci. Adv._ 2020, _6_ , eaax1950
  * Corso et al. 2023 Corso, G.; Stärk, H.; Jing, B.; Barzilay, R.; Jaakkola, T. DiffDock: Diffusion Steps, Twists, and Turns for Molecular Docking. International Conference on Learning Representations (ICLR). 2023
  * Han and Kim 2002 Han, K.-H.; Kim, J.-H. Quantum-inspired evolutionary algorithm for a class of combinatorial optimization. _IEEE Transactions on Evolutionary Computation_ 2002, _6_ , 580–593
  * Tang 2019 Tang, E. A Quantum-Inspired Classical Algorithm for Recommendation Systems. Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing. New York, NY, USA, 2019; p 217–228
  * Montiel Ross 2020 Montiel Ross, O. H. A Review of Quantum-Inspired Metaheuristics: Going From Classical Computers to Real Quantum Computers. _IEEE Access_ 2020, _8_ , 814–838
  * Arrazola et al. 2020 Arrazola, J. M.; Delgado, A.; Bardhan, B. R.; Lloyd, S. Quantum-inspired algorithms in practice. _Quantum_ 2020, _4_ , 307
  * Oshiyama and Ohzeki 2022 Oshiyama, H.; Ohzeki, M. Benchmark of quantum-inspired heuristic solvers for quadratic unconstrained binary optimization. _Sci. Rep._ 2022, _12_ , 2146
  * Kowalsky et al. 2022 Kowalsky, M.; Albash, T.; Hen, I.; Lidar, D. A. 3-regular three-XORSAT planted solutions benchmark of classical and quantum heuristic optimizers. _Quantum Science and Technology_ 2022, _7_ , 025008
  * Huang et al. 2022 Huang, T.; Xu, J.; Luo, T.; Gu, X.; Goh, R.; Wong, W.-F. Benchmarking Quantum(-inspired) Annealing Hardware on Practical Use Cases. 2022
  * Zeng et al. 2024 Zeng, Q.-G.; Cui, X.-P.; Liu, B.; Wang, Y.; Mosharev, P. A.; Yung, M.-H. Performance of Quantum Annealing Inspired Algorithms for Max-Cut Problem. 2024; in preparation
  * Cao et al. 2018 Cao, Y.; Romero, J.; Aspuru-Guzik, A. Potential of quantum computing for drug discovery. _IBM J. Res. Dev._ 2018, _62_ , 6:1–6:20
  * Outeiral et al. 2021 Outeiral, C.; Strahm, M.; Shi, J.; Morris, G. M.; Benjamin, S. C.; Deane, C. M. The prospects of quantum computing in computational molecular biology. _WIREs Comput. Mol. Sci._ 2021, _11_ , e1481
  * Zinner et al. 2021 Zinner, M.; Dahlhausen, F.; Boehme, P.; Ehlers, J.; Bieske, L.; Fehring, L. Quantum computing’s potential for drug discovery: Early stage industry dynamics. _Drug Discov. Today_ 2021, _26_ , 1680–1688
  * Blunt et al. 2022 Blunt, N. S.; Camps, J.; Crawford, O.; Izsák, R.; Leontica, S.; Mirani, A.; Moylett, A. E.; Scivier, S. A.; Sünderhauf, C.; Schopf, P.; Taylor, J. M.; Holzmann, N. Perspective on the Current State-of-the-Art of Quantum Computing for Drug Discovery Applications. _J. Chem. Theory Comput._ 2022, _18_ , 7001–7023
  * Mato et al. 2022 Mato, K.; Mengoni, R.; Ottaviani, D.; Palermo, G. Quantum molecular unfolding. _Quantum Sci. Technol._ 2022, _7_ , 035020
  * Pandey et al. 2022 Pandey, M.; Zaborniak, T.; Melo, H.; Galda, A.; Mulligan, V. K. Multibody molecular docking on a quantum annealer. 2022; arXiv:2210.11401
  * Yan et al. 2023 Yan, G.; Wu, H.; Yan, J. Quantum 3D Graph Learning with Applications to Molecule Embedding. Proceedings of the 40th International Conference on Machine Learning. 2023; pp 39126–39137
  * Robert et al. 2021 Robert, A.; Barkoutsos, P. K.; Woerner, S.; Tavernelli, I. Resource-efficient quantum algorithm for protein folding. _npj Quantum Information_ 2021, _7_ , 38
  * Babbush et al. 2014 Babbush, R.; Perdomo-Ortiz, A.; O’Gorman, B.; Macready, W.; Aspuru-Guzik, A. _Advances in Chemical Physics: Volume 155_ ; John Wiley & Sons, Ltd, 2014; pp 201–244
  * Perdomo-Ortiz et al. 2012 Perdomo-Ortiz, A.; Dickson, N.; Drew-Brook, M.; Rose, G.; Aspuru-Guzik, A. Finding low-energy conformations of lattice protein models by quantum annealing. _Sci. Rep._ 2012, _2_ , 571
  * Fingerhuth et al. 2018 Fingerhuth, M.; Babej, T.; Ing, C. A quantum alternating operator ansatz with hard and soft constraints for lattice protein folding. 2018; arXiv:1810.13411
  * Babej et al. 2018 Babej, T.; Ing, C.; Fingerhuth, M. Coarse-grained lattice protein folding on a quantum annealer. 2018; arxiv:1811.00713
  * Zhang 2011 Zhang, G. Quantum-inspired evolutionary algorithms: a survey and empirical study. _J. Heuristics_ 2011, _17_ , 303–351
  * Liu et al. 2023 Liu, B.; Wang, K.; Xiao, D.; Yu, Z. On connection among quantum-inspired algorithms of the Ising model. _Communications in Mathematical Sciences_ 2023, _21_ , 2013–2028
  * Goto 2016 Goto, H. Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network. _Sci. Rep._ 2016, _6_ , 21686
  * Goto et al. 2019 Goto, H.; Tatsumura, K.; Dixon, A. R. Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems. _Sci. Adv._ 2019, _5_ , eaav2372
  * Goto 2019 Goto, H. Quantum Computation Based on Quantum Adiabatic Bifurcations of Kerr-Nonlinear Parametric Oscillators. _J. Phys. Soc. Japan_ 2019, _88_ , 061015
  * Goto et al. 2021 Goto, H.; Endo, K.; Suzuki, M.; Sakai, Y.; Kanao, T.; Hamakawa, Y.; Hidaka, R.; Yamasaki, M.; Tatsumura, K. High-performance combinatorial optimization based on classical mechanics. _Sci. Adv._ 2021, _7_ , eabe7953
  * Kanao and Goto 2022 Kanao, T.; Goto, H. Simulated bifurcation for higher-order cost functions. _Appl. Phys. Express_ 2022, _16_ , 014501
  * Chermoshentsev et al. 2021 Chermoshentsev, D. A.; Malyshev, A. O.; Esencan, M.; Tiunov, E. S.; Mendoza, D.; Aspuru-Guzik, A.; Fedorov, A. K.; Lvovsky, A. I. Polynomial unconstrained binary optimisation inspired by optical simulation. 2021; arXiv:2106.13167
  * 49 See Supplemental Material for the details.
  * Gardner 1986 Gardner, M. _Knotted Doughnuts and Other Mathematical Entertainments_ ; Springer New York: New York: W. H. Freeman, 1986
  * Harris and Harris 2015 Harris, S.; Harris, D. _Digital Design and Computer Architecture: ARM Edition_ , 1st ed.; Morgan Kaufmann Publishers Inc.: San Francisco, CA, USA, 2015
  * 52 RDKit: Open-source cheminformatics. https://www.rdkit.org
  * Wang et al. 2023 Wang, B.-Y.; Cui, X.; Xu, H.; Zhang, G.; Yung, M.-H. Phase encoding method and device. 2023; Patent No.CN202310184040.5
  * Wójcikowski et al. 2015 Wójcikowski, M.; Zielenkiewicz, P.; Siedlecki, P. Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field. _J. Cheminformatics_ 2015, _7_ , 26

Figure 4: For Table of Contents Only


$f_{\mu}(\theta)=\sum_{\mathtt{s},\mathtt{r}}[R^{\mathrm{T}}\sigma_{\mu}R]_{\mathtt{s}\mathtt{r}}[\tilde{R}^{t-2}\tilde{\sigma}_{\mu}]_{\mathtt{s}\mathtt{r}}R_{\mathtt{r}\mathtt{s}}$
for any local pairing, where
$\tilde{\sigma}_{\mathtt{s}\mathtt{r}}=R_{\mathtt{s}\mathtt{r}}\sigma_{\mathtt{s}\mathtt{r}}$.
More explicitly, $f_{\mu}(\theta)$ for each $\mu$ reads

$f_{\mu}(\theta)=\begin{dcases*}\cos^{2}\theta(1+\cos^{t-2}2\theta)&$\mu=0$\\\
\sin^{2}\theta(1-\cos^{t-2}2\theta-2\cos^{t-1}2\theta)&$\mu=1$\\\
\sin^{2}\theta(1+\cos^{t-2}2\theta)&$\mu=2$\\\
\cos^{2}\theta(1-\cos^{t-2}2\theta+2\cos^{t-1}2\theta)&$\mu=3$\end{dcases*}.$
(95)

Note that $\sum_{\mu}f_{\mu}(\theta)=2(1+\cos^{t}2\theta)$, as expected.
Therefore any $s\neq 0$ pairing that contains a class I operator
$\sigma_{\mu_{x}}\otimes\tilde{P}_{\alpha_{x}}$ yields the factor
$q^{-1}f_{\mu_{x}}(\theta)$ at site $x$. This means that the contributions
from the spin and the colour degrees of freedom are again decoupled as in the
PSFF. The other two cases can be treated altogether as their actions on the
colour degree of freedom are trivial. First, we notice that there is only one
relevant local paring at $q\to\infty$ after Haar averaging, which is the one
associated with the identity permutation that carries the weight $q$. This
implies that the diagrams that include $\sigma_{1}$ or $\sigma_{2}$ have
vanishing contributions at large $q$. On the other hand, non-zero
contributions are produced when $\mu=0,3$ with the factor
$g_{\mu}=\sum_{\mathtt{s},\mathtt{r}}[R^{\mathrm{T}}\sigma_{\mu}R]_{\mathtt{s}\mathtt{s}}\tilde{R}^{t-1}_{\mathtt{s}\mathtt{r}}[\sigma_{\mu}]_{\mathtt{r}\mathtt{r}}$.
A straightforward calculation gives $g_{0}=2$ and $g_{3}=2\cos^{t}2\theta$.
Again, they sum to $2(1+\cos^{t}2\theta)$ as they should.

Now we are in the position to evaluate autocorrelation functions. Combining
the above results, we find that $\overline{C_{\alpha\alpha}(t)}$ behaves at
large $q$ as

$\displaystyle\overline{C_{\alpha\alpha}(t)}$
$\displaystyle=(2q)^{-L}q^{L-a}q^{-a}g_{0}^{L-b-a}g_{3}^{b}$
$\displaystyle\quad\times\prod_{\mu=0}^{3}f_{\mu}^{\mathtt{a}_{\mu}}(\theta)\prod_{m=1}^{n}\left[(TD_{\mathrm{c}})^{a_{m}}T\right]_{00}$
$\displaystyle=C^{\mathrm{s}}_{\alpha}(t)C^{\mathrm{c}}_{\alpha}(t)$ (96)

where
$C^{\mathrm{s}}_{\alpha}(t)=2^{-a}\cos^{bt}2\theta\prod_{\mu=0}^{3}f_{\mu}^{\textsection\mathtt{a}_{\mu}}(\theta)$
and $C^{\mathrm{c}}_{\alpha}(t)=q^{-2a}e^{-2n\varepsilon
t}(t-1)^{n}\left[1+(t-2)e^{-\varepsilon t}\right]^{a-n}$ Note that
$C^{\mathrm{c}}_{\alpha}(t)$ is nothing but $\overline{C_{\alpha\alpha}(t)}$
in the RPM without the spin degree of freedom Eq. (25). Since
$C^{\mathrm{s}}_{\alpha}(t)$ is a function that decays monotonically in time,
the contribution from the spin degree of freedom turns out to influence only
the physics after the peak, which comes from $C^{\mathrm{c}}_{\alpha}(t)$. To
be more precise, using the fact that asymptotically
$C^{\mathrm{s}}_{\alpha}(t)$ goes as $C^{\mathrm{s}}_{\alpha}(t)\simeq
2^{-a}\cos^{bt}2\theta$ provided that $\theta\neq 0$, the late time
asymptotics of the auto-correlation functions is given by

$(2q^{2})^{-a}\cos^{bt}2\theta e^{-2n\varepsilon t}(t-1)^{n}.$ (97)

This suggests that the operator strings that contribute to the decay of the
PSFF are i) those that are made of one cluster of class I operators and ii)
those that have only one cluster of class II operators with $b=1$.

Another important observation is that the decay time $t_{\mathrm{d}}$ becomes
nontrivial when the operator string is made of class II operators only, i.e.
$a=n=1$, in which case autocorrelation functions become
$\overline{C_{\alpha\alpha}(t)}=\cos^{bt}2\theta$. Invoking the sum rule Eq.
(5), it is thus these operators which give rise to the nontrivial value of the
PSFF at $t=1$ in Eq. (92).

Finally, let us also make a brief comment on the off-diagonal correlation
functions $\overline{\langle C_{\alpha\beta}(t)\rangle}$. Following the same
argument as above, it is readily seen that they vanish unless
$\tilde{P}_{\mu,\alpha}$ and $\tilde{P}_{\mu^{\prime},\beta}$ act on the
colour degrees of freedom in the same way, i.e. $P_{\alpha}=P_{\beta}$. The
nontrivial contribution comes from the Pauli matrices, and in particular, if
$\tilde{P}_{\alpha}$ and $\tilde{P}_{\beta}$ contain $\sigma_{\mu}$ and
$\sigma_{\mu^{\prime\prime}}$ on the same site, they amount to the factor
$f_{\mu\mu^{\prime}}=\sum_{\mathtt{s},\mathtt{r}}(R^{\mathrm{T}}\sigma_{\mu}R)_{\mathtt{s}\mathtt{r}}(\tilde{R}^{t-2}\tilde{\sigma}_{\mu^{\prime}})_{\mathtt{s}\mathtt{r}}R_{\mathtt{r}\mathtt{s}}$.
Note $f_{\mu\mu}=f_{\mu}$ and the only nonvanishing nondiagonal matrix
elements $f_{\mu\nu}$ are $f_{02}=-f_{20}=-{\rm i}\sin
2\theta(1-\cos^{t-2}2\theta)/2$ and $f_{13}=f_{31}=-\sin
2\theta(1+\cos^{t-2}2\theta+2\cos^{t-1})/2$. The net factor coming from the
Pauli matrix on each site can be obtained by multiplying $f_{\mu\nu}$.

###  A.3 OTOC of the RPM with the spin sector

As in the computation for the autocorrelation functions as well as the PSFF,
most of the techniques we developed in calculating the large-$q$ OTOC in Sec.
IV carry over even in the presence of the spin degree of freedom. We
parameterise the local operator as $\mathcal{O}=\sigma_{\mu}\otimes
P_{\alpha}$, where both $\sigma_{\mu}$ and $P_{\alpha}$ act on the same site.
Note that, as in the spinless case, the result does not depend on the choice
of $P_{\alpha}$ after Haar averaging.

The task again boils down to computing the on-site diagrams so let us start
with the sites that are neither at $0$ nor $x$.

####  A.3.1 Diagrams without local operators

The main change due to the additional on-site spin degree of freedom is that
now each on-site unitary is represented by a pair of black and white dots
together with a red square. While the dots act on the colour degree of freedom
as before, red squares correspond to the matrix $R$ that acts on the spin
degree of freedom (see Fig. 20a and 20b). With this rule it is straightforward
to see that, upon Haar averaging, the same leading pairings as in the OTOC of
the spinless RPM labeled by $p=0,\cdots,2t$ contribute, but in this case we
also need to account for a possible factor stemming from $R$ (see Fig. 20c-20f
for the possible local pairings at $t=1$. As in the spinless case, only Fig.
20c, 20d, and 20e contribute at large $q$). It turns out that this factor is
trivial due to the trace structure of the diagram, and it only gives the
factor $2$. Therefore, barring the sign, each contributing local pairing
carries the factor $2q$. The situation is, however, slightly subtler when the
diagram contains the operator $\mathcal{O}$, which we turn to next.

(a)  (b)

(c)  (d)  (e)

(f)

Figure 20: (a,b): On-site unitaries with the extra spin structure. (c,d,e,f):
Local pairings at $t=1$ with the factors (c,d) $2q$, (e): $-2q$, and (f):
$-2q^{-2}$ at large $q$. A red square indicates the matrix $R$.

####  A.3.2 Diagrams with local operators

In contrast to the case for sites on which no local operator acts, at sites
with local operators Haar averaging generates different local pairings
depending on whether the local operator $\mathcal{O}$ is class I or class II.
While in the former the situation is the same as in the spinless RPM, the case
where the operator is class II is a little more involved. Haar averaging in
fact gives a similar result to that of on-site diagrams on sites without local
operators. This is because the operator acts on the colour degree of freedom
trivially when the operator is class II. Hence we have $2t+1$ contributing
local pairings with the same transfer matrix $S$ for these diagrams as
previously. The difference, however, is that since the operator $\mathcal{O}$
still acts on the spin space nontrivially, every local pairing carries the
same factor $q\eta_{\mu}(\theta)$ where
$\eta_{\mu}(\theta)=\sum_{\mathtt{s}}[(R^{\mathrm{T}}\sigma_{\mu}R)_{\mathtt{s}\mathtt{s}}]^{2}$
except the one corresponding to $p=0$ pairing (on site $0$) and $p=2t$ pairing
(on site $x$). The factor $\eta_{\mu}(\theta)$ can be easily computed as
$\eta_{0}(\theta)=2,\eta_{1}(\theta)=2\sin^{2}\theta,\eta_{2}(\theta)=0,\eta_{3}(\theta)=2\cos^{2}\theta$.

####  A.3.3 Computation of the OTOC

Since there is essentially no effect of having the spin degree of freedom when
the local operator is a class I operator, the OTOC for a class I operator
$\mathcal{O}$ is exactly the same as the one we already computed Eq. (58).
When the operator is class II instead, the OTOC can be expressed as a trace

$\overline{\mathcal{C}(x,t)}=1-\frac{1}{4}\mathrm{Tr}(SG_{1}\hat{S}^{|x|-1}SG_{2}\hat{S}^{L-|x|-1}),$
(98)

where $G_{1}$ and $G_{2}$ are $(2t+1)\times(2t+1)$ diagonal matrices whose
diagonal entries are
$[G_{1}]_{pp}=2\delta_{0p}+(1-\delta_{0p})\eta_{\mu}(\theta)$ and
$[G_{2}]_{pp}=2\delta_{tp}+(1-\delta_{tp})\eta_{\mu}(\theta)$. It turns out
that we can evaluate (98) in a way that is similar to the case of class I
operators, which yields

$\overline{\mathcal{C}(x,t)}=\left(1-\frac{1}{2}\eta_{\mu}(\theta)\right)^{2}\sum_{m=0}^{t-|x|-1}\begin{pmatrix}t-1\\\
m\end{pmatrix}\rho^{m}(1-\rho)^{t-m-1},$ (99)

Thus at the coarse-grained scale (i.e. large $x$ and $t$), it behaves as

$\overline{\mathcal{C}(x,t)}=\left(1-\frac{1}{2}\eta_{\mu}(\theta)\right)^{2}\Phi\left(\frac{v_{B}t-|x|}{\sqrt{2Dt}}\right).$
(100)

Even though the OTOC now depends on $\theta$ explicitly, it is striking that
the way it does is rather simple and the additional spin space in the RPM only
induces the time-independent multiplicative factor fixed by
$\eta_{\mu}(\theta)$.

##  Appendix B Numerical methods

To evaluate numerically the various quantities considered in this work, we
rely on a matrix-free Monte Carlo method. First we describe the numerical
method used to determine the PSFF. The scheme requires storing and operating
on a complex vector of dimension $q^{L_{A}+L}=D^{2}_{A}D_{\bar{A}}$, and so is
less memory-intensive than working with the full Floquet operator for
$L_{A}<L$. As usual, we denote the two subregions in the system by $A$ and
$\bar{A}$, with $L_{\bar{A}}=L-L_{A}$. We also introduce a supplementary
subsystem $R$ consisting of $L_{A}$ sites. We initialise $A$ and $R$ in a
maximally entangled state $|\Phi_{AR}\rangle$ which we can view as a tensor
product of $L_{A}$ Bell pairs, each having one end in $A$ and the other end in
$R$. The region $\bar{A}$ is initialised in a Haar-random state
$|\psi\rangle_{\bar{A}}$. The full initial state is

$\displaystyle|\Phi_{AR}\rangle\otimes|\psi_{\bar{A}}\rangle=\Big{(}D_{A}^{-1/2}\sum_{i=0}^{D_{A}-1}|i\rangle_{R}\otimes|i\rangle_{A}\Big{)}\otimes|\psi\rangle_{\bar{A}}.$

We then evolve the system $A\bar{A}$ under the Floquet operator $W$. After $t$
time steps the state is

$\displaystyle 1_{R}\otimes
W^{t}\big{(}|\Phi_{AR}\rangle\otimes|\psi_{\bar{A}}\rangle\big{)}$ (101)
$\displaystyle=D_{A}^{-1/2}\sum_{i=0}^{D_{A}-1}|i\rangle_{R}\otimes
W^{t}\big{(}|i\rangle_{A}\otimes|\psi\rangle_{\bar{A}}\big{)}.$

Contracting this state with $\langle\Phi_{AR}|$ generates a pure state on
$\overline{A}$,

$\displaystyle D_{A}\langle\Phi_{AR}|1_{R}\otimes
W^{t}\big{(}|\Phi_{AR}\rangle\otimes|\psi_{\bar{A}}\rangle\big{)}=\big{[}\text{Tr}_{A}W^{t}\big{]}|\psi\rangle_{\bar{A}}.$

Finally we calculate the squared norm of the vector on the right-hand side.
Averaging this over $|\psi\rangle$ we generate

$\displaystyle\text{E}_{\psi}\Big{[}\langle\psi|_{\bar{A}}\big{[}\text{Tr}_{A}W^{-t}\big{]}\big{[}\text{Tr}_{A}W^{t}\big{]}|\psi\rangle_{\bar{A}}\Big{]}$
(102)
$\displaystyle=D_{\bar{A}}^{-1}\text{Tr}_{\bar{A}}\big{[}\text{Tr}_{A}W^{t}\text{Tr}_{A}W^{-t}\big{]}=K_{A}(t).$

Note that there are a number of different Monte Carlo methods that could be
designed to determine the PSFF; the advantage of this one is that we directly
contract the final-time indices in $\bar{A}$ between the forward and backward
evolution, and this drastically reduces the statistical fluctuations. The
variance over $|\psi\rangle$ is

$\displaystyle\text{Var}_{\psi}\Big{[}\langle\psi|_{\bar{A}}\big{[}\text{Tr}_{A}W^{-t}\big{]}\big{[}\text{Tr}_{A}W^{t}\big{]}|\psi\rangle_{\bar{A}}\Big{]}$
(103)
$\displaystyle=D_{\bar{A}}^{-2}\text{Tr}_{\bar{A}}\big{[}\text{Tr}_{A}W^{t}\text{Tr}_{A}W^{-t}\text{Tr}_{A}W^{t}\text{Tr}_{A}W^{-t}\big{]}.$

For $W$ a Haar-random $D\times D$ matrix, the average over the right-hand side
can in the large-$D$ limit be evaluated as $2D_{\bar{A}}^{-1}$. This is the
contribution from the Gaussian diagram with equal-time pairing. This implies
that for any given realization of the circuit, the fluctuations of
$\langle\psi|_{\bar{A}}\big{[}\text{Tr}_{A}W^{-t}\big{]}\big{[}\text{Tr}_{A}W^{t}\big{]}|\psi\rangle_{\bar{A}}$
around its average (i.e. the PSFF $K_{A}(t)$) are exponentially small in
$L_{\bar{A}}$.

The method used to determine the autocorrelation functions for operators with
support within $A$ is closely related. After creating the state in Eq. (101)
we construct the reduced density matrix for the region $AR$, and then evaluate
the expectation value of e.g. $P_{\alpha,R}\otimes P_{\alpha,A}$. Here we can
view $P_{\alpha,R}$ as the same operator as $P_{\alpha,A}$ but acting at the
initial time. Averaging the result of the calculation over
$|\psi_{\bar{A}}\rangle$ we find
$D_{\bar{A}}\text{Tr}[P_{\alpha}(t)P_{\alpha}(0)]=D_{A}C_{\alpha\alpha}(t)$.
Because the ensembles of circuits we consider enjoy a statistical invariance
under single-site unitary operations it is sufficient to calculate
autocorrelation functions for operators that are diagonal in the
comptutational basis. For this reason we need only compute the diagonal
elements of the density matrix for $AR$, and hence the memory required for the
final step of this calculation is only $q^{2L_{A}}=D_{A}^{2}$.

We now describe our scheme for evaluating $K^{B}_{A}(t)$ and $N^{B}_{A}(t)$.
To calculate $K^{B}_{A}(t)$ we first create two states
$|\Psi_{f}\rangle=|\Phi_{AR}\rangle\otimes|\psi_{\bar{A},f}\rangle$ and
$|\Psi_{b}\rangle=|\Phi_{AR}\rangle\otimes|\psi_{\bar{A},b}\rangle$, with
$|\psi_{\bar{A},f}\rangle$ and $|\psi_{\bar{A},b}\rangle$ independent random
states. Recall that $|\Phi_{AR}\rangle$ is a maximally entangled state of
$L_{A}$ qubits in $A$ and $L_{R}=L_{A}$ qubits in $R$. The states
$|\Psi_{f}\rangle$ and $|\Psi_{b}\rangle$ are evolved forward in time by
$\lceil t/2\rceil$, and backward in time by $\lfloor t/2\rfloor$,
respectively. At this step we have generated

$\displaystyle|\Psi_{f}(\lceil t/2\rceil)\rangle$ $\displaystyle=(1_{R}\otimes
W^{\lceil t/2\rceil})|\Psi_{f}\rangle.$ (104) $\displaystyle|\Psi_{b}(\lfloor
t/2\rfloor)\rangle$ $\displaystyle=(1_{R}\otimes W^{-\lfloor
t/2\rfloor})|\Psi_{b}\rangle.$

We now divide the $L$ sites of the system into subsystems $B$ and $\bar{B}$,
consisting of $L_{B}$ and $L_{\bar{B}}=L-L_{B}$ sites, and contract indices of
$|\Psi_{f}(\lceil t/2\rceil)\rangle$ and $\langle\Psi_{b}(\lfloor
t/2\rfloor)|$ in the regions $R$ and $B$. The result is a set of uncontracted
indices in $\bar{B}$. Computing the squared norm of the resulting vector, and
averaging over $|\psi_{\bar{A},f}\rangle$ and $|\psi_{\bar{A},b}\rangle$, we
find $(D_{A}/D_{B})K^{B}_{A}(t)$. In the modulus-squared object, we can
understand the above contraction of indices in $R$ as generating the sum over
all operators $\alpha\in A$ which defines the PSFF. Similarly, the contraction
of indices in $B$ generates a resolution of the identity in the operator space
for this region (corresponding to the sum over $\beta$ in Eq. (80)). The
contraction in $\bar{B}$, on the other hand, has the effect of eliminating
operator paths which, at time $\lceil t/2\rceil$, have nontrivial support in
$\bar{B}$.

To calculate $N^{B}_{A}(t)$, which in contrast with $K^{B}_{A}(t)$ is a sum
over probabilities of operator paths, we construct a pair of states
$|\Psi_{f}(\lceil t/2\rceil)\rangle$ and $|\Psi_{f^{\prime}}(\lceil
t/2\rceil)\rangle$, having identical structure but independent realisations of
the randomness (i.e. a statistically independent initial state of $\bar{A}$).
We then contract indices of $|\Psi_{f}(\lceil t/2\rceil)\rangle$ and
$\langle\Psi_{f^{\prime}}(\lceil t/2\rceil)|$ in the regions $R$ and $B$, and
compute the squared norm of the resulting state. After averaging over initial
states of $\bar{A}$, we obtain $N^{B}_{A}(t)$.

## References

  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, J. High Energy Phys. 2016, 1 (2016).
  * Nahum  _et al._ [2017] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys. Rev. X 7, 031016 (2017).
  * Nahum  _et al._ [2018] A. Nahum, S. Vijay, and J. Haah, Phys. Rev. X 8, 021014 (2018).
  * von Keyserlingk  _et al._ [2018] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, and S. L. Sondhi, Phys. Rev. X 8, 021013 (2018).
  * Fisher  _et al._ [2023] M. P. A. Fisher, V. Khemani, A. Nahum, and S. Vijay, Annu. Rev. Condens. Matter Phys. 14, 335 (2023).
  * Rakovszky  _et al._ [2018] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk, Phys. Rev. X 8, 031058 (2018).
  * Khemani  _et al._ [2018] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev. X 8, 031057 (2018).
  * Chan  _et al._ [2018a] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. X 8, 041019 (2018a).
  * Bertini  _et al._ [2018] B. Bertini, P. Kos, and T. Prosen, Phys. Rev. Lett. 121, 264101 (2018).
  * Mehta [2004] M. L. Mehta,  _Random matrices_ (Elsevier, 2004).
  * Joshi  _et al._ [2022] L. K. Joshi, A. Elben, A. Vikram, B. Vermersch, V. Galitski, and P. Zoller, Phys. Rev. X 12, 011018 (2022).
  * Chan  _et al._ [2018b] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. Lett. 121, 060601 (2018b).
  * Garratt and Chalker [2021a] S. J. Garratt and J. T. Chalker, Phys. Rev. Lett. 127, 026802 (2021a).
  * Prange [1997] R. E. Prange, Phys. Rev. Lett. 78, 2280 (1997).
  * Kunz [1999] H. Kunz, J. Phys. A: Math. Gen. 32, 2171 (1999).
  * Cotler  _et al._ [2017] J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, J. High Energy Phys. 2017, 1 (2017).
  * Gharibyan  _et al._ [2018] H. Gharibyan, M. Hanada, S. H. Shenker, and M. Tezuka, J. High Energy Phys. 2018, 124 (2018).
  * Braun  _et al._ [2020] P. Braun, D. Waltner, M. Akila, B. Gutkin, and T. Guhr, Phys. Rev. E 101, 052201 (2020).
  * Garratt and Chalker [2021b] S. J. Garratt and J. T. Chalker, Phys. Rev. X 11, 021051 (2021b).
  * Flack  _et al._ [2020] A. Flack, B. Bertini, and T. Prosen, Phys. Rev. Research 2, 043403 (2020).
  * Bertini  _et al._ [2021] B. Bertini, P. Kos, and T. Prosen, Commun. Math. Phys. 387, 597 (2021).
  * Claeys and Lamacraft [2020] P. W. Claeys and A. Lamacraft, Phys. Rev. Res. 2, 033032 (2020).
  * Rampp  _et al._ [2023] M. A. Rampp, R. Moessner, and P. W. Claeys, Phys. Rev. Lett. 130, 130402 (2023).
  * Thouless [1977] D. J. Thouless, Phys. Rev. Lett. 39, 1167 (1977).
  * Altshuler and Shklovskii [1986] B. L. Altshuler and B. I. Shklovskii, JETP 64, 127 (1986).
  * Bertrand and García-García [2016] C. L. Bertrand and A. M. García-García, Phys. Rev. B 94, 144201 (2016).
  * Sierant  _et al._ [2020] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev. Lett. 124, 186601 (2020).
  * Serbyn  _et al._ [2017] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. B 96, 104201 (2017).
  * Sonner  _et al._ [2021] M. Sonner, M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. B 104, L081112 (2021).
  * Nahum  _et al._ [2022] A. Nahum, S. Roy, S. Vijay, and T. Zhou, Phys. Rev. B 106, 224310 (2022).
  * Rakovszky  _et al._ [2022] T. Rakovszky, C. W. von Keyserlingk, and F. Pollmann, Phys. Rev. B 105, 075131 (2022).
  * von Keyserlingk  _et al._ [2022] C. von Keyserlingk, F. Pollmann, and T. Rakovszky, Phys. Rev. B 105, 245101 (2022).
  * Swingle  _et al._ [2016] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, Phys. Rev. A 94, 040302 (2016).
  * Vermersch  _et al._ [2019] B. Vermersch, A. Elben, L. M. Sieberer, N. Y. Yao, and P. Zoller, Phys. Rev. X 9, 021061 (2019).
  * Lieb and Robinson [1972] E. H. Lieb and D. W. Robinson, Commun. Math. Phys. 28, 251 (1972).
  * McCulloch and von Keyserlingk [2022] E. McCulloch and C. W. von Keyserlingk, J. Phys. A: Math. Theor. 55, 274007 (2022).
  * Brouwer and Beenakker [1996] P. W. Brouwer and C. W. J. Beenakker, J. Math. Phys. 37, 4904 (1996).
  * Akila  _et al._ [2016] M. Akila, D. Waltner, B. Gutkin, and T. Guhr, J. Phys. A: Math. Theor. 49, 375101 (2016).
  * Lerose  _et al._ [2021] A. Lerose, M. Sonner, and D. A. Abanin, Phys. Rev. X 11, 021040 (2021).
  * Friedman  _et al._ [2019] A. J. Friedman, A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. Lett. 123, 210603 (2019).
  * Roy and Prosen [2020] D. Roy and T. Prosen, Phys. Rev. E 102, 060202 (2020).
  * Moudgalya  _et al._ [2021] S. Moudgalya, A. Prem, D. A. Huse, and A. Chan, Phys. Rev. Research 3, 023176 (2021).
  * Parker  _et al._ [2019] D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, Phys. Rev. X 9, 041017 (2019).
  * Chan  _et al._ [2021] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. Research 3, 023118 (2021).
  * Bertini  _et al._ [2019] B. Bertini, P. Kos, and T. Prosen, Phys. Rev. Lett. 123, 210601 (2019).
  * Gutkin  _et al._ [2020] B. Gutkin, P. Braun, M. Akila, D. Waltner, and T. Guhr, Phys. Rev. B 102, 174307 (2020).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Gopalakrishnan  _et al._ [2015] S. Gopalakrishnan, M. Müller, V. Khemani, M. Knap, E. Demler, and D. A. Huse, Phys. Rev. B 92, 104202 (2015).
  * Crowley and Chandran [2022] P. J. D. Crowley and A. Chandran, SciPost Phys. 12, 201 (2022).
  * Garratt  _et al._ [2021] S. J. Garratt, S. Roy, and J. T. Chalker, Phys. Rev. B 104, 184203 (2021).
  * Garratt and Roy [2022] S. J. Garratt and S. Roy, Phys. Rev. B 106, 054309 (2022).
  * Long  _et al._ [2023] D. M. Long, P. J. D. Crowley, V. Khemani, and A. Chandran, Phys. Rev. Lett. 131, 106301 (2023).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. Van Den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, _et al._ , Nature 618, 500 (2023).
  * Kechedzhi  _et al._ [2023] K. Kechedzhi, S. Isakov, S. Mandrà, B. Villalonga, X. Mi, S. Boixo, and V. Smelyanskiy, Future Gener. Comput. Syst. , (2023).
  * Kos  _et al._ [2018] P. Kos, M. Ljubotina, and T. Prosen, Phys. Rev. X 8, 021062 (2018).
  * Liao and Galitski [2022] Y. Liao and V. Galitski, Phys. Rev. Research 4, L012037 (2022).