compat=1.1.0

# Information scrambling and butterfly velocity in quantum spin glass chains

Venkata Lokesh K. Y1, Surajit Bera2 and Sumilan Banerjee2 1Centre of High
Energy Physics, Indian Institute of Science, Bangalore 560012, India
2Centre for Condensed Matter Theory, Department of Physics, Indian Institute
of Science, Bangalore 560012, India<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1Centre of High Energy Physics, Indian
Institute of Science, Bangalore 560012, India
2Centre for Condensed Matter Theory, Department of Physics, Indian Institute
of Science, Bangalore 560012, India

###### Abstract

We make lattice generalization of two well-known zero dimensional models of
quantum spin glass, Sachdev-Ye (SY) and spherical quantum $p$-spin glass
model, to one dimension for studying crossovers in non-local scrambling
dynamics due to glass transition, complex dynamics, and quantum and thermal
fluctuations in paramagnetic (PM) and spin glass (SG) phases. In the SY chain
of quantum dots, each described by infinite-range random Heisenberg model of
$N$ spin-$S$ $SU(M)$ spins, we obtain the quantum Lyapunov exponent
$\lambda_{\mathrm{L}}$ and butterfly velocity $v_{B}$ as a function of
temperature $T$ and the quantum parameter $S$ across the PM-SG phase boundary
using a bosonic spinon representation in the large-$N,M$ limit. In particular,
we extract asymptotic $T$ and $S$ dependence, e.g., power laws, for
$\lambda_{\mathrm{L}}$ and $v_{B}$ in different regions deep inside the phases
and near the replica symmetry breaking SG transition. We find the chaos to be
non-maximal almost over the entire phase diagram. Very similar results for
chaos indicators are found for the $p$-spin glass chain as a function of
temperature and a suitable quantum parameter $\Gamma$, with some important
qualitative differences. In particular, $\lambda_{\mathrm{L}}$ and $v_{B}$
exhibit a maximum, coinciding with onset of complex glassy relaxation, above
the glass transition as a function of $T$ and $\Gamma$ in the PM phase of the
$p$-spin glass model. In contrast, the maximum is only observed as a function
of $S$, but not with temperature, in the PM phase of SY model. The maximum
originates from enhanced chaos due to maximal complexity in the glassy
landscape. Thus, the results in the SY model indicate very different evolution
of glassy complexity with quantum and thermal fluctuations.

##  1 Introduction

In a generic quantum many-body system, quantum information encoded in local
degrees of freedom in the initial state _scrambles_ to spread globally over
the entire system due to its time evolution. The dynamics of quantum
information, e.g., how rapidly local information scrambles, characterized by a
rate or the quantum Lyapunov exponent $\lambda_{\mathrm{L}}$, and how fast the
information propagates from one spatial region to other, quantified by a speed
or the butterfly velocity $v_{B}$, can provide crucial insights into the non-
trivial dynamics of plethora of quantum many-body systems Sekino and Susskind
(2008); Maldacena  _et al._ (2016); Kitaev ; Sahu  _et al._ (2019). These
insights are complementary to that from transport and usual dynamical
correlations. As a result, information scrambling has become an important
theme of research for understanding holographic dualityJensen (2016), black
hole information scramblingSekino and Susskind (2008); Maldacena  _et al._
(2016), non-Fermi liquids Kitaev ; Maldacena and Stanford (2016); Patel and
Sachdev (2017); Banerjee and Altman (2017); Chowdhury  _et al._ (2022) and
quantum thermalization Sahu  _et al._ (2019); Srednicki (1994); Murthy and
Srednicki (2019). For systems close to a semiclassical limit, information
scrambling can often be connected to chaos Aleiner  _et al._ (2016); Stanford
(2016); Rozenbaum  _et al._ (2017); Scaffidi and Altman (2019); Xu  _et al._
(2020), namely extreme sensitivity of the classical trajectories to initial
conditions.

For systems with many metastable configurations, information scrambling via
quantum dynamics may become enhanced when the dynamics samples these multitude
of configurations, leading to extreme sensitivity to initial conditions in the
classical limit. Glasses Kob (1997); Debenedetti and Stillinger (2001);
Berthier and Biroli (2011), which possess extensive number of metastable
states Debenedetti and Stillinger (2001); Berthier and Biroli (2011); Büchner
and Heuer (1999); Sastry (2001); Sciortino (2005), are among the most
prominent examples of such systems. At low temperatures, the system remains
trapped near one of the metastable configurations for long times and explore
the phase space extremely slowly. As a result, glasses exhibit many
fascinating dynamical properties ranging from complex relaxation with
anomalously large time scales Binder and Young (1986); Kob (1997), aging and
memory effects Vincent (2007), unusual transport properties such as the
violation of Stokes-Einstein relations Hodgdon and Stillinger (1993),
dynamical heterogeneity Berthier and Biroli (2011); Karmakar  _et al._ (2010)
etc. In particular, spin glasses Edwards and Anderson (1975); Binder and Young
(1986); Mezard  _et al._ (1987) with quenched randomness undergo ergodicity
breaking thermodynamic phase transition at low temperature, e.g., for mean-
field or infinite-range models, to replica symmetry broken Parisi (1979);
Castellani and Cavagna (2005) spin glass phase. Recent work Bera  _et al._
(2022) on a well-known solvable zero-dimensional or infinite-range model of
quantum spin glass, namely the quantum $p$-spin glass model Nieuwenhuizen
(1995); Nieuwenhuizen and Ritort (1998); Cugliandolo  _et al._ (2000), have
shown that chaos, as characterized by the Lyapunov exponent
$\lambda_{\mathrm{L}}$, becomes stronger over a region above the glass
transition where the dynamics exhibits typical signatures Kob (1997) of
complex glassy relaxation. In this region of enhanced chaos, the dynamics
samples exponentially large number of saddle points of the classical potential
energy landscape and becomes maximally complex Correale  _et al._ (2023);
Auffinger and Ben Arous (2011), before eventually getting trapped within the
basin of attraction of one of the glassy minima below the glass transition. In
this work, we study scrambling dynamics in two solvable, albeit _spatially
extended models_ , of quantum spin glasses and show how the measures of
information scrambling, $v_{B}$ and $\lambda_{\mathrm{L}}$, are influenced by
glass transition, complex dynamics, and quantum and thermal fluctuations in
various paramagnetic and spin glass phases.

Figure 1.1: An one-dimensional (1D) generalization of (a) SY model, (b)
$p$-spin model for $p=3$. Each lattice site contains a quantum dot described
by zero-dimensional SY model in (a) and quantum spherical $p$-spin glass model
in (b) with $N\gg 1$ spins (blue circles). In (a), the coupling between
nearest neighbour sites are random two-spin interaction with one from each
lattice site for the SY model. In the case of $p$-spin model in (b), the
random nearest neighbour interaction is between two spins from one site and
one form the other and vice versa.

In line with the growing interests Bera  _et al._ (2022); Anous and Haehl
(2021); Correale  _et al._ (2023) in understanding the relation between
scrambling and the complex dynamics of spin glass (SG) systems, we consider
two well-known large-$N$ zero dimensional ($0+1$ D) quantum spin glass models,
namely (a) Sachdev-Ye (SY) model, a random Heisenberg spin glass Sachdev and
Ye (1993); Georges  _et al._ (2000, 2001), and (b) quantum spherical $p$-spin
glass model Nieuwenhuizen (1995); Nieuwenhuizen and Ritort (1998); Cugliandolo
_et al._ (2000), and make lattice generalizations of these models to one
dimension ($1+1$ D), as shown schematically in Fig.1.1. These generalizations
allow us to meaningfully define a butterfly velocity $v_{B}$, along with
$\lambda_{\mathrm{L}}$, from an out-of-time-order correlation (OTOC) function
Kitaev ; Gu  _et al._ (2017). The SY model is parent to the celebrated
fermionic Sachdev-Ye-Kitaev (SYK) model Kitaev ; Maldacena and Stanford
(2016); Sachdev (2015), which has emerged as a paradigm to study quantum many-
body chaos, information scrambling Gu  _et al._ (2017); Banerjee and Altman
(2017); Chowdhury  _et al._ (2022), as well as strongly correlated metals Song
_et al._ (2017a). The SYK model is maximally chaotic, namely it saturates the
remarkable upper bound Maldacena  _et al._ (2016); Kitaev ; Maldacena and
Stanford (2016) $2\pi k_{\mathrm{B}}T/\hbar$ for $\lambda_{\mathrm{L}}$. There
have been numerous lattice and higher dimensional generalizations Gu  _et al._
(2017); Song  _et al._ (2017b); Davison  _et al._ (2017); Jian  _et al._
(2017); Jian and Yao (2017); Berkooz  _et al._ (2017); Zhang (2017); Chowdhury
_et al._ (2018); Haldar  _et al._ (2018) of the zero dimensional SYK model,
with a view of studying transport and information propagation in extended
interacting systems. Our lattice generalizations of the spin glass models pave
the way to access non-local dynamical correlations and spatio-temporal
information scrambling in spatially extended strongly disordered systems that
undergo ergodicity-breaking glass transitions.

Earlier works on quantum chaos in spin glass systemsBera  _et al._ (2022);
Anous and Haehl (2021); Correale  _et al._ (2023) have only looked at the
Lyapunov exponent, e.g., its dependence on temperature $T$ and quantum
fluctuation parameter in the paramagnetic and spin glass phases, in the zero
dimensional quantum $p$-spin glass model, i.e., for a single _$p$ -spin glass
quantum dot_ (see Fig.1.1). On the other hand, studies on spatially extended
large-$N$ systems, such as lattice generalizations of the SYK model Gu  _et
al._ (2017); Song  _et al._ (2017b); Davison  _et al._ (2017); Jian  _et al._
(2017); Jian and Yao (2017); Berkooz  _et al._ (2017); Zhang (2017); Chowdhury
_et al._ (2018); Haldar  _et al._ (2018), have been limited to non (replica)
symmetry broken phases, such as Fermi and non-Fermi liquids metals. Our
lattice generalizations of the spin glass models retain the same large-$N$
saddle points as the zero dimensional models, and capture replica symmetry
broken spin glass (SG) and unbroken paramagnetic (PM) phases, such as a spin
liquid for SY quantum dot. Thus, the spatial locality in these extended models
permits us to explore information scrambling across replica-symmetry breaking
spin glass phase transition. In particular, we obtain the variations of the
Lyapunov exponent and butterfly velocity with thermal and quantum fluctuations
across the PM-SG phase transition and various crossovers between dynamically
different paramagnetic states, e.g., spin liquid and local moment states in
the SY model. We also find the asymptotic dependence of $\lambda_{\mathrm{L}}$
and $v_{B}$ on temperature $T$ and quantum/classical parameter in different
regions deep inside the phases as well as near the replica symmetry breaking
SG transitions.

Temperature dependence of chaos indicators, especially that of the butterfly
velocity, has been studied in several earlier works in various contexts, for
examples, for solvable large-$N$ models Gu  _et al._ (2017); Guo  _et al._
(2019); Chowdhury and Swingle (2017); Sahu and Swingle (2020), perturbatively
in interaction for weakly interacting systems Patel  _et al._ (2017) and in
$1/N$ expansion in strongly interacting systems Patel and Sachdev (2017), via
numerical methods Sahu and Swingle (2020) as well as in the semiclassical
limit of spin models Bilitewski  _et al._ (2018); Ruidas and Banerjee (2021);
Bilitewski  _et al._ (2021). For a maximally chaotic
[$\lambda_{\mathrm{L}}=2\pi T$ ($\hbar,k_{\mathrm{B}}=1$)] interacting
diffusive metal built from SYK dots Gu  _et al._ (2017), $v_{B}$ has been
found to follow a $\sqrt{T}$ dependence. Also, in this case, the quantity,
$v_{B}^{2}/\lambda_{\mathrm{L}}$, which has the dimension of the diffusion
constant, turns out to be exactly equal to energy diffusion constant, as
expected from holographic theories Blake (2016). Similar relation,
$D=v_{B}^{2}/4\lambda_{L}$ has been found to be satisfied for electron
diffusion constant $D$ in weakly interacting diffusive metal Swingle and
Chowdhury (2017); Patel  _et al._ (2017), where $\lambda_{\mathrm{L}}\sim
T^{d/2},~{}v_{B}\sim T^{d/4}$ in $d$ dimension, and a temperature independent
$D$ is determined by elastic scattering. A power-law $T$ dependence of $v_{B}$
is similar to that at a critical point Wei  _et al._ (2019); Ruidas and
Banerjee (2021), i.e., $v_{B}\sim T^{1-1/z}$ for a dynamical exponent $z>1$.
In a strongly interacting system with critical Fermi surface of $N$-species
fermions coupled with $U(1)$ gauge field, power-law behaviours,
$\lambda_{\mathrm{L}}\sim T$, $v_{B}\sim T^{1/3}$, with energy diffusion
constant $\sim v_{B}^{2}/\lambda_{\mathrm{L}}$, have been obtained Patel and
Sachdev (2017) at leading order in $1/N$-expansion. The Lyapunov exponent and
butterfly velocity also have been computed in (2+1 d) $O(N)$ non-linear sigma
model, where, in the quantum critical region Chowdhury and Swingle (2017),
$\lambda_{\mathrm{L}}\sim T/N$, and $v_{B}\sim c$, a constant close to the
speed of light in the theory. Similar $T$ independent $v_{B}$ is found in the
$O(N)$-symmetry broken ordered phase Chowdhury and Swingle (2017) with
$\lambda_{\mathrm{L}}\sim T^{3}$, whereas $\lambda_{\mathrm{L}}\sim
e^{-E_{g}/T}$, due to the gap $E_{g}$ in the disordered phase for $T\ll
E_{g}$, and $v_{B}$ obeys a $\sqrt{T}$ behaviour Chowdhury and Swingle (2017);
Sahu and Swingle (2020). The latter is consistent with numerical calculations
in 1D quantum Ising model Sahu and Swingle (2020). Power-law temperature
dependence of $\lambda_{\mathrm{L}}$ and $v_{B}$ is also observed Bilitewski
_et al._ (2018); Ruidas and Banerjee (2021) in the ordered and disordered
phases of spin models with semiclassical dynamics.

Figure 1.2: Schematic phase diagram for chaos in SY model. Dependence of
Lyapunov exponent $\lambda_{\mathrm{L}}$ and butterfly velocity $v_{B}$ on
temperature $T$ and quantum parameter $S$ are shown in the paramagnetic (PM)
and spin glass (SG) phases. $\lambda_{\mathrm{L}}$ and $v_{B}$ increases as a
power law in $S$ for small $S$ in the PM phase, and as a power law in $T$ for
low temperature in the SG phase. Both the chaos indicators increases slowly,
as $\sim\ln{T}$, above the glass transition $T_{SG}$, before saturating to a
constant at high temperature. In the SG phase, $\lambda_{\mathrm{L}}$
decreases as $\sim 1/S^{2}$, whereas $v_{B}\sim e^{aS}$ increases
exponentially in $S$. There is also maximum (not shown) in both
$\lambda_{\mathrm{L}}$ and $v_{B}$ as a function of $S$ in the PM phase, at
intermediate $S$, but there are no maxima as a function of $T$ (see main
text). The region bounded by $T_{SG}(S)$ and $T/\tilde{J}=1$ ($\tilde{J}$
coupling constant) line (dotted line) demarcates the quantum critical region
within the PM phase. Arrows indicate the direction of increase of the
parameter. Figure 1.3: Schematic phase diagram for chaos in $p$-spin glass
model. Dependence of $v_{B}$ on temperature $T$ and quantum parameter $\Gamma$
are shown in the PM and SG phases. Similar dependence has been seen for
$\lambda_{\mathrm{L}}$ in Ref.Bera  _et al._ , 2022. The continuous PM-SG
phase transition $T_{d}(\Gamma)$ [or $T_{SG}(\Gamma)$] (solid line) terminates
at tricritical point and thereafter the transition is first order (hatched
line). There is a maximum in $v_{B}$ at $T_{m}(\Gamma)$ (dotted line), both as
a function of $T$ and $\Gamma$. $v_{B}$ increases as power law with $T$ and
$\Gamma$ at low temperature and small $\Gamma$, respectively, in the SG phase.
The butterfly velocity increases as $\sim\ln{T}$ and as a power law in
$\Gamma$ above the glass transition in the PM phase, before the maximum.
Arrows indicate the direction of increase of the parameter.

Our results for the dependence of $\lambda_{\mathrm{L}}$ and $v_{B}$ on
temperature $T$ and the quantum/classical parameter $S$ in the chain of SY
quantum dots are summarized on the schematic phase diagram of Fig.1.2. In this
model, the dots are described by infinite-range random Heisenberg model of
spin-$S$ $SU(M)$ spins on $N$ sites, where $N,M\to\infty$. Like in the lattice
generalization Gu  _et al._ (2017) of the SYK model, the SY chain model, in a
bosonic spinon representation, is described by the same saddle point of the
zero-dimensional quantum dots with the strength of the coupling $J$
renormalized by the inter-dot interaction. As a result,
$\lambda_{\mathrm{L}}$, which is extracted from a local OTOC, has the same $S$
and $T$ dependence of a single dot, or the original SY model Sachdev and Ye
(1993). On the contrary, the butterfly velocity can only be defined in a
lattice model like the SY chain of Fig.1.1(a). To compare with the chaotic
properties of the SY chain, we also calculate the butterfly velocity from the
non-local OTOC as a function of $T$ and quantum parameter $\Gamma$ (defined
later) in a similar lattice model [Fig.1.1(b)] of quantum spherical $p$-spin
spin glass model Nieuwenhuizen (1995); Nieuwenhuizen and Ritort (1998);
Cugliandolo  _et al._ (2000), as summarized in Fig.1.3 (schematic). Again,
$\lambda_{\mathrm{L}}(\Gamma,T)$ of the $p$-spin glass chain is identical to
that of the corresponding zero-dimensional model, and was studied in ref.Bera
_et al._ , 2022. Our results are summarized below.

(1) As shown in Fig.1.2, we find a power-law increase with $S$ in the SY
model, i.e.,
$\lambda_{\mathrm{L}}(S,T)\simeq\lambda_{\mathrm{L}}(0,T)+aS^{\alpha_{\lambda}}$
and $v_{B}(S,T)\simeq v_{B}(0,T)+aS^{\alpha_{v}}$ for small $S$ in the PM
phase with $\alpha_{\lambda(v)}\approx 0.5$ at a fixed temperature, like in
the quantum critical region for $T<J$, where two types of PM behaviours, spin
liquid and local moment, compete. In contrast to the power-law or gap-like
temperature dependence seen in the disordered phase of many other models
Chowdhury and Swingle (2017); Sahu and Swingle (2020); Bilitewski  _et al._
(2018); Ruidas and Banerjee (2021), a $\ln{T}$ variation of both
$\lambda_{\mathrm{L}}$ and $v_{B}$ are seen in the PM phase pf the SY model
close to the transition temperature $T_{SG}$. Similar $\ln{T}$ dependence is
seen for the PM phase of the $p$-spin glass model, where
$v_{B}\sim\Gamma^{\alpha_{v}}$ with $\alpha_{v}\approx 1$.

(2) On the contrary, in the SG phase of SY model, power-law temperature
dependence, $\lambda_{\mathrm{L}}(S,T)\sim T^{\beta_{\lambda}}$, and a power-
law increment $v_{B}(S,T)\simeq v_{B}(S,0)+aT^{\beta_{v}}$, with
$\beta_{\lambda}\approx 1-1.5$ and $\beta_{v}\approx 1.5$, are observed. The
velocity $v_{B}(S,0)$ extrapolated to zero temperature is finite, indicating
information spreading in the quantum SG ground state of the SY chain. Very
similar behaviours of $\lambda_{\mathrm{L}}(\Gamma,T)$ and $v_{B}(\Gamma,T)$
occurs in the SG phase of $p$-spin glass. In the SG phase of SY model at a
fixed $T$, the Lyapunov exponent and butterfly velocity exhibit contrasting
trends with increasing $S$ or decreasing quantum fluctuations;
$\lambda_{\mathrm{S}}$ decreases as $\sim 1/S^{2}$, whereas $v_{B}$ increases
exponentially with $S$. These behaviours are different from the $p$-spin glass
model where both $\lambda_{\mathrm{L}}$ Bera  _et al._ (2022) and $v_{B}$
decreases with decreasing quantum fluctuation, e.g., $v_{B}(\Gamma,T)\simeq
v_{B}(0,T)+a\Gamma^{\alpha_{v}}$ with $\alpha_{v}\approx 0.5$.

(3) The Lyapunov exponent and butterfly velocity have smooth crossovers across
the PM-SG transitions in both SY and $p$-spin glass models. However, $v_{B}$
exhibits a dip as a function of $S$ at the SG transition for the SY chain,
indicating slowing down of information spread at the transition. This is
unlike the $v_{B}(\Gamma)$ which monotonically changes across the SG
transition in the $p$-spin glass.

(4) As mentioned earlier, $\lambda_{\mathrm{L}}(\Gamma,T)$ exhibits a maximum
Bera  _et al._ (2022) [$T_{m}(\Gamma)$ in Fig.1.3] in the PM phase of the
$p$-spin glass model as function of both $\Gamma$ and $T$. The maximum
corresponds to enhanced chaos due to the sampling of maximum complexity of the
glassy landscape at $T_{m}$ Correale  _et al._ (2023); Auffinger and Ben Arous
(2011). The same maximum is found in $v_{B}$ for $p$-spin glass. However, the
maximum in $\lambda_{\mathrm{L}}$ and $v_{B}$ is seen in the PM phase of the
SY model only as a function of $S$, but not as a function of $T$. This implies
very different evolution of complexity with quantum and thermal fluctuations
in SY model, unlike that in the $p$-spin glass model.

(5) For computing the butterfly velocity, we implement a numerical approach
alternative to the more standard semi-analytical _single-mode ansatz_ Gu and
Kitaev (2019); Guo  _et al._ (2019). The results from the numerical method
agree very well with those from the single-mode ansatz over the entire phase
diagrams of SY and $p$-spin glass chains, thus, independently verifying the
applicability of the single-mode ansatz even for the replica symmetry broken
spin glass phase. Using the single-mode ansatz, we find that chaos is non-
maximal almost over the entire phase diagram for SY and $p$-spin glass models.

We also connect our results for the chaos parameters with spectral
characteristics of various regions of the phase diagram of the two models,
like spin liquid, local moment, classical and quantum paramagnets, and spin
glass. The rest of the paper is organized as follows. In Sec.2, we discuss the
one dimensional lattice generalizations of the zero-dimensional SY and quantum
$p$-spin glass model. The large-$N$ saddle-point equations for the PM and SG
phased for both the models are described in Sec.3, followed by review of the
phase diagrams known from earlier works. Sec.4 describes the formalism for
calculating local OTOC and the Lyapunov exponent using ladder approximation
for the PM and SG phases for the SY model. This section reports the results
for $\lambda_{\mathrm{L}}(S,T)$ in the SY model and compare them with those in
the $p$-spin glass model. The formalism for computing the butterfly velocity
from non-local OTOC using two different methods, a numerical method and the
semi-analytical single-mode ansatz Gu and Kitaev (2019), is discussed in Sec.5
for both SY and $p$-spin glass models. The results for $v_{B}$ as a function
of $T$ and quantum parameter for the two models are described in this section.
The appendices (Appendix A, Appendix B) give additional details of the
derivations of the saddle-point equations and methods for their numerical
self-consistent solutions. Appendix A.2 discusses the results for the spectral
properties for the SY model. Appendix C provides details of the numerical
computation of $v_{B}$.

##  2 Lattice generalizations of solvable quantum spin glass models

As mentioned earlier, we consider one-dimensional (1D) chains of zero-
dimensional random quantum Heisenberg and spherical $p$-spin glass quantum
dots, as shown in Fig.1.1 and discussed below.

###  2.1 1D chain of Sachdev-Ye (SY) spin-glass quantum dots

We make a lattice generalization of the well-known zero-dimensional SY model
of random quantum Heisenberg model by arranging $x=1,\dots,L$ quantum dots of
SY model, where each dot contains $i=1,\dots,N$ sites. The model is described
by the Hamiltonian

$\displaystyle\mathcal{H}=\frac{1}{2\sqrt{MN}}\sum_{x,i\neq
j}J_{ij,x}S_{i\alpha\beta x}S_{j\beta\alpha x}$
$\displaystyle+\frac{1}{\sqrt{MN}}\sum_{x,i,j}J^{\prime}_{ij,x}S_{i\alpha\beta
x}S_{j\beta\alpha,x+1}.$ ( 2.1)

Here the first term represents all-to-all infinite-range interactions within a
dot, and the second term the nearest-neighbour inter-dot coupling. The
couplings $J_{ij,x}$ and $J^{\prime}_{ij,x}$ are drawn from Gaussian
distribution, independently at each site with zero mean and variances
$\overline{J_{ij}^{2}}=J^{2}$ and $\overline{J_{ij}^{\prime 2}}=J^{\prime 2}$.
$S_{i\alpha\beta x}$ denote the spin degrees of freedom, the generators of the
SU$(M)$ group for spin $S$, at each $i$ and dot $x$, where
$\alpha,\beta=1,\dots,M$. The model is exactly solvable in the large $N$ limit
($N\to\infty$) followed by the large $M$ limit ($M\to\infty$), i.e.,
thermodynamics, as well as equilibrium and non-equilibrium properties can be
exactly obtained in this limit through large-$N$ saddle-point equations, as
known from previous works Sachdev and Ye (1993); Georges  _et al._ (2000,
2001); Biroli and Parcollet (2002) on the zero-dimensional model. The effect
of quantum fluctuations in the model can be tuned by varying the spin $S$ and
the classical and semiclassical regimes can be probed by taking the large $S$
($S\to\infty$) limit while appropriately scaling the energy /temperature
scales Georges  _et al._ (2001). The $M=2$ limit of the model is the usual
$SU(2)$ Heisenberg model, which has been studied via quantum Mote Carlo (QMC)
Grempel and Rozenberg (1998) after taking the $N\to\infty$ limit. The large
$M$ limit was found to capture the properties of the $SU(2)$ model quite well.

In Sec.3.1, we discuss the equilibrium phase diagram obtained from the
large-$N$ saddle-point equations as a function of temperature $T$ and spin
$S$. The large-$N$ saddle point is obtained via the standard bosonic
representation of the $SU(M)$ spin Sachdev and Ye (1993). Subsequently, our
main results on the dynamical phase diagram for chaos, characterized via the
Lyapunov exponent $\lambda_{\mathrm{L}}$ and butterfly velocity $v_{B}$, are
discussed in Secs.4,5.

###  2.2 1D chain of spherical $p$-spin glass quantum dots

To compare and contrast the results for $\lambda_{\mathrm{L}}$ and $v_{B}$ in
the SY chain, we consider the lattice generalization of another paradigmatic
model of quantum glasses, the quantum spherical $p$-spin glass model
Nieuwenhuizen (1995); Nieuwenhuizen and Ritort (1998); Cugliandolo  _et al._
(2000, 2001). The 1D chain of the $p$-spin glass quantum dots ($x=1,\dots,L$)
is described by the Hamiltonian

$\displaystyle\mathcal{H}=\sum_{x,i}\frac{\pi_{i,x}^{2}}{2\mathcal{M}}+\frac{1}{3!}\sum_{x,ijk}J_{ijk,x}s_{i,x}s_{j,x}s_{k,x}$
$\displaystyle+\frac{1}{2!}\sum_{x,ijk}s_{i,x}(J^{+}_{ijk,x}s_{j,x+1}s_{k,x+1}+J^{-}_{ijk,x}s_{j,x-1}s_{k,x-1}),$
( 2.2)

for $p=3$. The intra- and inter-dot all-to-all couplings among sites
$i=1,\dots,N$ in the dots, $J_{ijk,x}$ and $J_{ijk,x}^{\pm}$, are Gaussian
random numbers with zero mean and variances $3J^{2}/(2N^{2})$ and $3J^{\prime
2}/(2N^{2})$, respectively. The commutation relation of the spin $s_{i,x}$
with the momentum $\pi_{j,x}$,
$[s_{i,x},\pi_{j,x^{\prime}}]=i\hbar\delta_{ij}\delta_{xx^{\prime}}$ induces
the quantum dynamics. The spherical constraint $\sum_{i}s^{2}_{i}=N$ makes the
model non trivial. Essentially, when the quantum dots are uncoupled
($J^{\prime}=0$), the Hamiltonian describes quantum particles with mass
$\mathcal{M}$ moving on the surface of an $N$-dimensional hypersphere. The
inter-dot interaction couples particles in neighbouring dots. The advantage of
this model lies in the fact that one can take a classical limit of the model
by continuously tuning a quantum fluctuation parameter
$\Gamma=\hbar^{2}/\mathcal{M}J$ to zero Bera  _et al._ (2022); Anous and Haehl
(2021). The classical model for $\hbar\to 0$ limit (fixed $\mathcal{M}$) for a
single dot gives rise dynamics identical to the mode coupling theory (MCT)
dynamics in the super-cooled liquid regime of structural glasses Götze (2008);
Kob (1997); Reichman and Charbonneau (2005); Berthier and Biroli (2011); Bera
_et al._ (2022). The limit $\mathcal{M}\to\infty$ while keeping $\hbar$ fixed
lead to a different classical limit of infinitely heavy mass Bera  _et al._
(2022); Anous and Haehl (2021). In this work we vary the dimensionless quantum
parameter $\Gamma$.

In Sec.3.2, we discuss the equilibrium phase diagram of the model as a
function of $T$ and $\Gamma$. The results for $\lambda_{\mathrm{L}}$ and
$v_{B}$ are discussed in Secs.4 and 5 and compared with those from the lattice
SY model.

##  3 Equilibrium phase diagrams of the spin glass chains

We first discuss the equilibrium phase diagram of the quantum spin glass
chains [Eq.(2.1) and Eq.(2.2)] as a function of $T$ and quantum fluctuation
parameters, $S$ or $\Gamma$ [Figs.3.1, 3.2]. The equilibrium phase diagrams
are obtained by solving large-$N$ saddle-point equations. As we discuss below,
due to the particular choice of inter-dot couplings in Eqs.(2.1),(2.2), the
large-$N$ saddle-point equations of the lattice models are identical to the
zero-dimensional saddle-point equations Sachdev and Ye (1993); Cugliandolo
_et al._ (2001) of the single dot with renormalized intra-dot coupling.

Figure 3.1: Mean-field phase diagram of the bosonic spinon model Georges  _et
al._ (2001) showing the paramagnetic and spin glass phases. Spin glass phase
below $T_{SG}$ (solid line) is determined by the _marginality_ criterion (see
main text). $T_{eq}$ (thick dashed line) is determined by equilibrium
criterion. The regimes within spin glass phase separated by blue lines are
identified via a large $S$ expansion of saddle point equations Georges  _et
al._ (2001) $I$, $II$ and $III$ are classical, semi-classical and quantum spin
glass regions respectively. $IV$ is the quantum critical regime in PM phase.
We perform the calculations of Lyapunov exponent and butterfly velocity across
the SG-PM transition for several $S$ as a function of $T$ (along vertical
dashed lines $a$, $b$ and $c$) and for several $T$ as a function of $S$ (along
horizontal dashed lines $d$, $e$ and $f$).

###  3.1 Sachdev-Ye chain

To capture both (replica) symmetry broken ordered SG and disordered
paramagnetic (PM) phases we use bosonic spinons to represent spin-$S$ SU$(M)$
generators as $S_{i\alpha\beta x}=b_{i\alpha x}^{\dagger}b_{i\beta
x}-S\delta_{\alpha\beta}$, where $b_{i\alpha x}$, $b_{i\alpha x}^{\dagger}$
($\alpha=1,\dots,M$) are bosonic annihilation and creation operators with the
constraint $\sum_{\alpha}b_{i\alpha x}^{\dagger}b_{i\alpha x}=SM$ ($S>0$) on
the number of bosons at each $i,x$ to fix the spin of the representation to
$S$. The saddle-point equations, their solutions and resulting phase diagram
for the bosonic representation have been discussed in detail for the zero-
dimensional model in the pioneering works Sachdev and Ye (1993); Georges  _et
al._ (2000, 2001). We briefly discuss the derivation of saddle point equations
for the SY chain model [Eq.(2.1)] in Appendix A.

The disorder averaged Green’s function
$G_{x}^{ab}(\tau,\tau^{\prime})=G^{ab}(\tau,\tau^{\prime})=-\frac{1}{M}\sum_{\alpha}\langle
T_{\tau}b_{i\alpha x}^{a}(\tau)\bar{b}_{i\alpha x}^{b}(\tau^{\prime})\rangle$
is independent of $i,x$ and determines the equilibrium dynamics of the model.
The Green’s function is obtained from the saddle point equations in the limit
$N\rightarrow\infty$ followed by $M\rightarrow\infty$ as

$\displaystyle(G^{-1})^{ab}(\mathfrak{i}\omega_{k})$
$\displaystyle=\left(\mathfrak{i}\omega_{k}+\lambda\right)\delta_{ab}-\Sigma^{ab}(\mathfrak{i}\omega_{k})$
( 3.1a) $\displaystyle\Sigma^{ab}(\tau)$
$\displaystyle=\tilde{J}^{2}[G^{ab}(\tau)]^{2}G^{ba}(-\tau)$ ( 3.1b)

where $\omega_{k}=2\pi kT$ is bosonic Matsubara frequency with $k$ an integer
($k_{B}=1$), $T$ the temperature in units of $\tilde{J}$. We have introduced
the replicas $a=1,\dots,n$ to carry out the disorder average.
$G(\mathfrak{i}\omega_{k})=\int_{0}^{\beta}d\tau
e^{\mathfrak{i}\omega_{k}\tau}G(\tau)$ and $\beta$ is the inverse temperature.
Lagrange multiplier, $\lambda$ is introduced to impose the constraint,
$G^{aa}(\tau=0^{-})=-S$. The above saddle-point equations are same as the
zero-dimensional ones Sachdev and Ye (1993); Georges  _et al._ (2000) with
renormalized coupling $\tilde{J}^{2}=J^{2}+J^{\prime 2}$ due to the nearest-
neighbour inter-dot coupling in Eq.(2.1). In the limit $n\rightarrow 0$, the
replica structure of the Green’s function determines the PM and SG phases as
we discuss below.

####  3.1.1 Paramagnetic phase

In the PM phase one employs a replica diagonal and symmetric ansatz,
$G^{ab}(\tau)=G(\tau)\delta_{ab}$. Thus, the saddle-point equations reduce to,

$\displaystyle[G(\mathfrak{i}\omega_{k})]^{-1}$
$\displaystyle=(\mathfrak{i}\omega_{k}+\lambda)-\Sigma(\mathfrak{i}\omega_{k})$
( 3.2a) $\displaystyle\Sigma(\tau)$
$\displaystyle=\tilde{J}^{2}[G(\tau)]^{2}G(-\tau)$ ( 3.2b)

Analytic continuation from Matsubara to real frequency
$\mathfrak{i}\omega_{k}\rightarrow\omega+\mathfrak{i}0^{+}$ gives the saddle
point equation for retarded Green’s function

$\displaystyle[G_{R}(\omega)]^{-1}=\omega+\lambda-\Sigma_{R}(\omega)$ ( 3.3)

where
$\Sigma_{R}(\omega)=\Sigma(\mathfrak{i}\omega_{k}\to\omega+\mathfrak{i}0^{+})$
and $\Sigma(\mathfrak{i}\omega_{k})=\int_{0}^{\beta}d\tau
e^{\mathfrak{i}\omega_{k}\tau}\Sigma(\tau)$. $\lambda$ is determined by the
constraint on $G(\tau=0^{-})$. The procedure to numerically solve for the
retarded Green’s function is discussed in Appendix A.

####  3.1.2 Spin glass phase

Following earlier works Georges  _et al._ (2001), we describe the spin glass
phase via the following one-step replica symmetry breaking (1RSB) ansatz for
the bosonic Green’s function that appears in Eq.(3.1).

$\displaystyle G^{ab}(\tau)$
$\displaystyle=\tilde{G}(\tau)\delta_{ab}-g\epsilon_{ab}$ ( 3.4)

or,
$G^{ab}(\mathfrak{i}\omega_{k})=\tilde{G}(\mathfrak{i}\omega_{k})\delta_{ab}-\beta
g\epsilon_{ab}\delta_{\omega_{k},0}$, and

$\Sigma^{ab}(\tau)=\tilde{\Sigma}(\tau)\delta_{ab}-\tilde{J}^{2}g^{3}\epsilon_{ab},$
( 3.5)

where $\epsilon_{ab}=1$ in the diagonal block of size $m\times m$, and 0
otherwise, for the $n\times n$ matrices. Here $\tilde{G}(\tau)$ is the regular
part of the Green’s function such that $\tilde{G}(\tau\to\infty)\to 0$ at
$T=0$. The Edward-Anderson (EA) spin-glass order parameter is obtained from
$q_{\mathrm{EA}}=\lim_{\tau\to\infty}(1/M^{2})\sum_{\alpha\beta}\langle
S_{i\alpha\beta x}(\tau)S_{i\beta\alpha
x}(0)\rangle$$=\lim_{\tau\to\infty}G^{aa}(\tau)G^{aa}(-\tau)=g^{2}$ at zero
temperature. Using
$[G^{-1}]^{ab}(\mathfrak{i}\omega_{k})=A(\mathfrak{i}\omega_{k})\delta_{ab}+B(\mathfrak{i}\omega_{k})\epsilon_{ab}$
and the standard replica matrix inversion, we get

$\displaystyle A(\mathfrak{i}\omega_{k})$
$\displaystyle=\frac{1}{\tilde{G}(\mathfrak{i}\omega_{k})},$ ( 3.6a)
$\displaystyle B(\mathfrak{i}\omega_{k})$ $\displaystyle=\frac{\beta
g}{[\tilde{G}(\mathfrak{i}\omega_{k})-m\beta
g]\tilde{G}(\mathfrak{i}\omega_{k})}.$ ( 3.6b)

The saddle point equation for $\tilde{G}$ is therefore,

$\displaystyle[\tilde{G}(\mathfrak{i}\omega_{k})]^{-1}=\mathfrak{i}\omega_{k}-\frac{\tilde{J}g}{\Theta}-[\tilde{\Sigma}(\mathfrak{i}\omega_{k})-\tilde{\Sigma}(\mathfrak{i}\omega_{k}=0)]$
( 3.7a) where $\displaystyle\tilde{\Sigma}(\tau)=$
$\displaystyle\tilde{J}^{2}\Big{[}\tilde{G}^{2}(\tau)\tilde{G}(-\tau)-2g\tilde{G}(\tau)\tilde{G}(-\tau)-g\tilde{G}^{2}(\tau)$
$\displaystyle+2g^{2}\tilde{G}(\tau)+g^{2}\tilde{G}(-\tau)\Big{]}$ ( 3.7b)

In the above, we have eliminated $\lambda$ using the equation for
$\mathfrak{i}\omega_{k}=0$ and the parameterization
$\tilde{G}(\mathfrak{i}\omega_{k}=0)=-\Theta/\tilde{J}g$. The saddle point
equation for $g$ is obtained from Eq.(3.6b) as,

$m\beta=\frac{1}{\tilde{J}g^{2}}\Big{(}\frac{1}{\Theta}-\Theta\Big{)}$ ( 3.8)

The equations Eq.3.7a, Eq.3.7 and Eq.3.8 along with the constraint,
$\tilde{G}(\tau=0^{-})=g-S$ provide a closed set of self-consistent saddle-
point equations for the SG phase. The characteristic of 1RSB symmetry breaking
solution is that the saddle point equations form a one-parameter family,
parameterized here by $\Theta$ or equivalently by the breaking point $m$. To
obtain this parameter we impose the usual marginality criterion Georges  _et
al._ (2001), namely we fix $\Theta=\Theta_{R}$, and hence $m$, such that the
lowest eigenvalue of the fluctuation matrix in the replica space vanishes.

The marginality criterion leads to a gapless
$\chi^{\prime\prime}(\omega)\propto\omega$ with a delta function peak at
$\omega=0$, where
$\chi^{aa}(\tau)=G^{aa}(\tau)G^{aa}(-\tau)$$=(1/M^{2})\sum_{\alpha\beta}\langle
S_{i\alpha\beta x}(\tau)S_{i\beta\alpha x}(0)\rangle$ is the local spin
susceptibility Georges  _et al._ (2001). The equilibrium criterion, which
minimizes free energy with respect to the breaking point $m$, leads to a
gapped spectrum Georges  _et al._ (2001). Analytic continuation of Eq.(3.7a)
gives the equation for retarded Green’s function
$G_{R}^{ab}(\omega)=G_{R}(\omega)\delta_{ab}$ with

$\displaystyle[G_{R}(\omega)]^{-1}=\omega-\frac{\tilde{J}g}{\Theta_{R}}-[\Sigma_{R}(\omega)-\Sigma_{R}(\omega=0)]$
( 3.9)

where
$\Sigma_{R}(\omega)=\tilde{\Sigma}(\mathfrak{i}\omega_{k}\to\omega+\mathfrak{i}0^{+})$
and $\tilde{\Sigma}(\mathfrak{i}\omega_{k})=\int_{0}^{\beta}d\tau
e^{\mathfrak{i}\omega_{k}\tau}\tilde{\Sigma}(\tau)$. $\lambda$ is fixed by the
constraint $\tilde{G}(\tau=0^{-})=g-S$, $g$ is calculated from the relation
Eq.3.8. The procedure for numerically calculating the retarded Green’s
function $G_{R}(\omega)$ is discussed in Appendix A. For a fixed $S$, the
temperature $T_{SG}$, below which the spins freeze into a SG order is
identified by $m(T_{SG})=1$, where the SG order parameter $q_{\mathrm{EA}}$
jumps discontinuously to zero in the PM phase Georges  _et al._ (2001). A
similar procedure is followed to obtain $S_{SG}$ for a given temperature $T$.
Unlike earlier works Georges  _et al._ (2001); Camjayi and Rozenberg (2003) ,
in this work we mostly use the real-frequency saddle-point solutions to find
the Lagrange multiplier $\lambda$, EA order parameter $q_{\mathrm{EA}}$ in the
SG phase, the PM-SG phase boundary, as well as the retarded Green’s function
$G_{R}(\omega)$, and chaos diagnostics $\lambda_{\mathrm{L}}$ and
$v_{\mathrm{B}}$. The numerical methods employed for obtaining the real-
frequency saddle-point solutions are discussed in Appendix A.1.1. The
imaginary-time saddle-point equations are only used for deriving the real-
frequency saddle-point equations via analytical continuation, and
occasionally, to check the consistency of the thermodynamic quantities like
$\lambda,~{}q_{EA}$ and $T_{SG}(S)$ computed from the real-frequency solution.
For the sake of completeness, we discuss the numerical method used for the
imaginary-time solution in Appendix A.1.2. The same procedure is followed for
the $p$-spin glass model.

####  3.1.3 Phase transition and crossovers

The equilibrium phase diagram of the SY chain model is identical to the zero-
dimensional case, which has been studied in detail in the $N,M\to\infty$ limit
in both bosonic and fermionic representation of the $SU(M)$ spins for general
$S$ Sachdev and Ye (1993); Georges  _et al._ (2000, 2001). The phase diagram
has been also studied for finite $M$ in the $N\to\infty$ limit, using $1/M$
expansion Christos  _et al._ (2022) and for small $S$, e.g., $S=1/2$, through
QMC Grempel and Rozenberg (1998) in the $N\to\infty$ limit, as well as via
exact diagonalization for finite $N$ Arrachea and Rozenberg (2002). In
Fig.3.1, we reproduce the schematic large-$N,M$ phase diagram obtained via
bosonic representation in ref.Georges  _et al._ , 2001. This is consistent
with QMC result Grempel and Rozenberg (1998) for $S=1/2$ for $N\to\infty$. The
system undergoes a PM to SG phase transition at a temperature
$T_{SG}\sim(2/3\sqrt{3})\tilde{J}S^{2}$, where the transition temperature,
$T_{SG}^{eq}$ obtained from the equilibrium criterion is slightly smaller than
the temperature calculated via the marginal stability criterion. We refer to
this latter temperature as $T_{SG}$. This corresponds to a dynamical
transition, similar to that in the $p$-spin glass model Cugliandolo  _et al._
(2001); Nieuwenhuizen and Ritort (1998), where the spin relaxation time
diverges and glassy aging behaviour commences in the non-equilibrium dynamics
Biroli and Parcollet (2002), approaching from the PM side.

The point $S=0,~{}T=0$ corresponds to a quantum critical point with a spin
liquid ground state, akin to SYK non-Fermi liquid Kitaev ; Maldacena and
Stanford (2016). The spin liquid state exhibits so-called ‘marginal spectrum’
Sachdev and Ye (1993); Georges  _et al._ (2000, 2001), i.e.,
$\chi^{\prime\prime}(\omega)\sim\mathrm{sgn}(\omega)$ at $T=0$, or
$\chi^{\prime\prime}(\omega,T)\sim\tanh(\omega/2T)$ at finite temperature, and
diverging spin susceptibility $\chi(T)\sim\ln(T/\tilde{J})$, whereas the
spectral density of the bosons
$\rho(\omega)=-(1/\pi)\mathrm{Im}G_{R}(\omega)\sim-\mathrm{sgn}(\omega)/\sqrt{|\omega|}$
has a characteristic power-law divergence at $T=0$. At any finite $S$ the
ground state is a spin glass and the spin liquid state becomes metastable up
to a small maximum $S=S_{max}\simeq 0.052$ Georges  _et al._ (2001) at $T=0$.
However, as shown in Fig.3.1, a quantum critical regime with marginal spectrum
persists at finite temperature, $T_{SG}<T\
{\raise-2.15277pt\hbox{$\buildrel<\over{\sim}$}}\ \tilde{J}$, up to $S\
{\raise-2.15277pt\hbox{$\buildrel<\over{\sim}$}}\ 2$. Apart from the spin
liquid, another paramagnetic solution to the Eqs.(3.2), the so-called ‘local
moment’ solution Georges  _et al._ (2001) coexists in the phase diagram. This
leads to a high-temperature ($T\gg\tilde{J}$) Curie susceptibility
$\chi(T)\simeq S(S+1)/T$ which gradually decreases to $\chi(T)\simeq S^{2}/T$
as $T\to 0$, i.e., below $T_{SG}$, due to partial screening of the spins due
to spin-spin interaction. Of course, for $T<T_{SG}$, the local moment solution
is metastable. In appendix Appendix A.2, we discuss the numerically obtained
bosonic spectral function $\rho(\omega)$ across the PM-SG transition as a
function of $S$ for fixed $T$, and as a function of $T$ for fixed $S$, in the
regime of interest for our chaos calculations (Secs.Appendix A.1.2, 5).

Based on the large-$S$ expansions in ref.Georges  _et al._ , 2001, the SG
phase in Fig.3.1 can be divided into three regimes – $I$. a classical regime
for $\tilde{J}S\lesssim T\lesssim\tilde{J}S^{2}$ where all non-zero Matsubara
frequency components of the Green’s function in eq.(3.1) can be neglected and
spins become commuting vectors of length $S$, $II$. a semiclassical regime for
$\tilde{J}\sqrt{S}\lesssim T\lesssim\tilde{J}S$, where
$\rho(\omega)\sim\omega$ for $\omega\to 0$ is obtained by $1/S$ expansion of
the saddle-point equations in the SG phase (Sec.3.1.2) after scaling $\omega$
and $T$ by $\tilde{J}S$, and $III$. a quantum regime for $0\leq
T\lesssim\tilde{J}\sqrt{S}$, deduced from the $1/S$ expansion of the internal
energy and linear in $T$ specific heat. In this work, we do not study the
chaos across these various quantum-classical crossovers in detail. As shown in
Fig.3.1, we focus on changes in the chaotic properties, like
$\lambda_{\mathrm{L}}$ and $v_{B}$, from different regions of the PM phase to
the SG phase across the SG transition.

Here, we only study the chaos in the SY chain model [Eq.(2.1)] using bosonic
representation. The model can also be studied in the large $N,M$ limit via
fermionic representation Sachdev and Ye (1993); Georges  _et al._ (2000);
Christos  _et al._ (2022). In this case, the ground state is a spin liquid for
all $0\leq S\leq 1$ accessible in the fermionic representation. Thus system is
maximally chaotic with $\lambda_{\mathrm{L}}=2\pi T$ ($\hbar=1$) for $T\to 0$.
However, for finite $M$ the SG order can be recovered below
$T_{SG}\sim\tilde{J}\exp(-\sqrt{M\pi})$ Christos  _et al._ (2022); Georges
_et al._ (2001). In the SG ground state at $T=0$, one obtains the marginal
spectrum $\chi^{\prime\prime}(\omega)\sim\mathrm{sgn}(\omega)$ above a
characteristic frequency $\omega^{*}=\tilde{J}q_{\mathrm{EA}}$, i.e., over a
range $\omega^{*}<|\omega|\ll\tilde{J}$, and the marginal SG behaviour for
$\chi^{\prime\prime}(\omega)$ for smaller frequencies $0<|\omega|<\omega^{*}$.
The EA order parameter in the fermionic case is $q_{\mathrm{EA}}\sim 1/M$.

###  3.2 $p$-spin glass chain

As in the case of SY model, the equilibrium dynamics of the $p$-spin glass
chain [Eq.(2.2)] is described by the disorder averaged Green’s function
$Q_{ab}(\tau,\tau^{\prime})=(1/N)\sum_{i}\overline{\langle\mathcal{T}_{\tau}s_{i}^{a}(\tau)s_{i}^{b}(\tau^{\prime})\rangle}$.
This model is solvable in the large $N$ limit and is described by the Dyson
equation,

$\displaystyle(Q^{-1})^{ab}(\mathfrak{i}\omega_{k})$
$\displaystyle=\left(\omega_{k}^{2}/\Gamma+z\right)\delta_{ab}-\Sigma^{ab}(\mathfrak{i}\omega_{k})$
( 3.10a) $\displaystyle\Sigma^{ab}(\tau)$
$\displaystyle=\frac{3\tilde{J}^{2}}{2}[Q_{ab}(\tau)]^{2},$ ( 3.10b)

where $\tilde{J}^{2}=J^{2}+J^{\prime 2}$. The derivation of the above saddle-
point equations, which are identical to that of the zero-dimensional $p$-spin
glass model Cugliandolo  _et al._ (2001), is discussed in Appendix B. The
spherical constraint $Q_{aa}(\tau=0)=1$ is imposed by the Lagrange multiplier
$z$. Similar to the case of SY model, in the limit $n\rightarrow 0$, we have
the replica diagonal PM and 1RSB SG phases. In the PM phase, the saddle point
equations, after analytical continuation, simplify to

$\displaystyle[Q_{R}(\omega)]^{-1}=-\omega^{2}/\Gamma+z-\Sigma_{R}(\omega),$
$\displaystyle\Sigma_{R}(\omega)=\Sigma_{R}(\mathfrak{i}\omega_{k}\to\omega+\mathfrak{i}0^{+})$
( 3.11a)

In the SG phase, we consider the 1RSB solution, which is exact Cugliandolo
_et al._ (2001). Using this ansatz, we can write the Dyson equation as,

$\displaystyle
Q^{ab}(\tau)=(q_{d}(\tau)-q_{\mathrm{EA}})\delta_{ab}+q_{\mathrm{EA}}\epsilon_{ab}$
( 3.12)

which after Fourier transformation becomes,

$\displaystyle
Q^{ab}(\mathfrak{i}\omega_{k})=(q_{d}(\mathfrak{i}\omega_{k})-\tilde{q}_{EA})\delta_{ab}+\tilde{q}_{EA}\epsilon_{ab},$
( 3.13)

where $\tilde{q}_{\mathrm{EA}}=\beta q_{\mathrm{EA}}\delta_{\omega_{k},0}$.
Further, it is convenient to write
$q_{d}(\tau)=\tilde{Q}(\tau)+q_{\mathrm{EA}}$, i.e., in terms of a regular
part, such that $\tilde{Q}(\tau\to\infty)\to 0$ at $T=0$, and EA order
parameter $q_{\mathrm{EA}}$. This can be inverted in the replica space as,

$\displaystyle[Q^{-1}]^{ab}(\mathfrak{i}\omega_{k})=A(\mathfrak{i}\omega_{k})\delta_{ab}+B(\mathfrak{i}\omega_{k})\epsilon_{ab}$
( 3.14)

with

$\displaystyle
A(\mathfrak{i}\omega_{k})=\frac{1}{q_{d}(\mathfrak{i}\omega_{k})-\tilde{q}_{\mathrm{EA}}},$
( 3.15a) $\displaystyle
B(\mathfrak{i}\omega_{k})=\frac{-\tilde{q}_{\mathrm{EA}}}{q_{d}(\mathfrak{i}\omega_{k})^{2}+(m-2)\tilde{q}_{EA}q_{d}(\mathfrak{i}\omega_{k})-(m-1)\tilde{q}_{\mathrm{EA}}^{2}}$
( 3.15b)

where $q_{\mathrm{EA}}$ is Edward-Anderson order parameter and
$\epsilon_{ab}=1$ for the diagonal blocks and zero otherwise.

This leads to the saddle point equations,

$\displaystyle\omega_{k}^{2}/\Gamma+z=\frac{1}{\tilde{Q}(\mathfrak{i}\omega_{k})}+\tilde{\Sigma}(\mathfrak{i}\omega_{k}),$
( 3.16)

and

$\displaystyle\frac{3(\beta
m)^{2}}{2}q_{\mathrm{EA}}^{3}=\frac{x_{p}^{2}}{1+x_{p}},$ ( 3.17a)
$\displaystyle y={\beta q_{\mathrm{EA}}}/{q_{d}(0)}\quad\text{and}\quad
x_{p}={my}/{(1-y)}.$ ( 3.17b)

where
$\tilde{\Sigma}(\mathfrak{i}\omega_{k})=\frac{3\tilde{J}^{2}}{2}\int_{0}^{\beta}d\tau
e^{\mathfrak{i}\omega_{k}\tau}(q_{d}^{2}(\tau)-q_{\mathrm{EA}}^{2})$. The
marginality criterion leads to $x_{p}=1$ Cugliandolo  _et al._ (2001), which
determines the breaking point, $m$. The saddle-point equations along with the
constraint $\tilde{Q}(\tau=0)=1-q_{\mathrm{EA}}$ , can be analytically
continued to real frequency and solved numerically. The details are discussed
in the Appendix B. We show the phase diagram determined by the solutions of
the saddle point equation in Fig.3.2.

Figure 3.2: Mean-field phase diagram of quantum spherical $p$-spin glass model
Cugliandolo  _et al._ (2001) as a function of temperature $T$ and
dimensionless quantum fluctuation parameter $\Gamma$ (see main text). The
solid line denotes second order phase transition between SG and PM phases. The
second order line terminates at a tricritical point and the transition becomes
first order for low temperature and larger $\Gamma$. Spin glass phase is
determined by the marginality criterion. We calculate the Lyapunov exponent
and butterfly velocity across the SG-PM transition for two values of $T$ as a
function of $\Gamma$ (along vertical dashed lines $x$ and $y$) and for two
values of $\Gamma$ as a function of $T$ (along horizontal dashed lines $u$ and
$v$).

##  4 Out-of-time-ordered correlator (OTOC) and Lyapunov exponent

To characterize the many-body quantum chaos in the SY and $p$-pin glass chain
models [Eqs.(2.1),(2.2)], we first extract Lyapunov exponent from _on-site_
OTOC. The calculation of Lyapunov exponent in the chain models are identical,
with the effective coupling $\tilde{J}=\sqrt{J^{2}+J^{\prime 2}}$, to that in
the corresponding zero-dimensional models, as have been done already for the
$p$-spin glass in Ref.Bera  _et al._ , 2022. Thus, we focus on the effective
zero-dimensional SY model, where the dot indices $x$ in Eq.(2.1) drop out, and
compare the results for $\lambda_{\mathrm{L}}(T)$ across the SG transition
with that in the zero-dimensional $p$-spin glass model Bera  _et al._ (2022).

Figure 4.1: Schwinger-Keldysh contour with four real-time branches separated
by $\beta/4$ in imaginary time, used for calculating OTOC.

To diagnose the quantum chaos in the SY model, we consider two disorder-
averaged four-point ‘regularized’ OTOCs Kitaev ; Maldacena and Stanford
(2016); Aleiner  _et al._ (2016) for the bosonic operators on a Schwinger
Keldysh (SK) contour Fig.4.1, with two forward and two backward branches,
equally spaced in imaginary time with separation $\beta/4$, namely

$\displaystyle F_{1}(t_{1},t_{2})$
$\displaystyle=\frac{1}{N^{2}M^{2}}\sum_{ij,\alpha\beta}\overline{\mathrm{Tr}[yb^{\dagger}_{i\alpha}(t_{1})yb^{\dagger}_{j\beta}(0)yb_{i\alpha}(t_{2})yb_{j\beta}(0)]}$
( 4.1a) $\displaystyle F_{2}(t_{1},t_{2})$
$\displaystyle=\frac{1}{N^{2}M^{2}}\sum_{ij,\alpha\beta}\overline{\mathrm{Tr}[yb_{i\alpha}(t_{1})yb^{\dagger}_{j\beta}(0)yb_{i\alpha}^{\dagger}(t_{2})yb_{j\beta}(0)]},$
( 4.1b)

where
$y^{4}=\mathrm{exp}(-\beta\mathcal{H})/\mathrm{Tr}[\mathrm{exp}(-\beta\mathcal{H})]$.
In this work, we do not consider the OTOC of the physical spin operators
$S_{i\alpha\beta x}$. Calculations of such an OTOC would involve 8-point
bosonic correlation functions. We assume the maximum Lyapunov exponent in the
large-$N,M$ SY model to be the same for OTOCs for bosonic and spin operators.

The calculation of the OTOCs in Eq.(4.1) on the SK contour (Fig.4.1) can be
done in the replica diagonal PM phase following a procedure very similar to
SYK and related models Kitaev ; Maldacena  _et al._ (2016); Banerjee and
Altman (2017). However below the dynamical transition temperature $T_{SG}$, in
the SG phase, calculating any dynamic correlation function is more involved
since the system always stays out of equilibrium Cugliandolo and Lozano
(1999); Biroli and Parcollet (2002). Alternatively, as discussed in Ref.Bera
_et al._ , 2022 for $p$-spin glass, the OTOCs can be calculated in the
marginal SG phase by replicated generating function Bera  _et al._ (2022)
$\mathcal{Z}^{n}$ on the four-branch contour [Fig.4.1], i.e.,

$\displaystyle\mathcal{Z}^{n}=$
$\displaystyle\frac{1}{Z^{n}}\int\mathcal{D}(\bar{b}_{i\alpha}^{a},b_{i\alpha}^{a})\mathcal{D}\lambda_{i}^{a}\mathrm{exp}(\mathfrak{i}S)$
( 4.2)

where $a=1,\dots,n$ are the replica indices and the action

$\displaystyle\mathfrak{i}S=$
$\displaystyle\Bigg{[}\mathfrak{i}\int_{\mathcal{C}}dz\sum_{\mathfrak{i}\alpha
a}\left\\{\bar{b}_{i\alpha}^{a}\mathfrak{i}\partial_{z}b_{i\alpha}^{a}-\lambda_{i}^{a}\left(\bar{b}_{i\alpha}^{a}(z)b_{i\alpha}^{a}(z)-SM\right)\right\\}$
$\displaystyle-\frac{1}{\sqrt{MN}}\frac{1}{2!}\mathfrak{i}\int_{\mathcal{C}}dz\sum_{ij,\alpha\beta,a}J_{ij}\bar{b}_{i\alpha}^{a}(z)b_{i\beta}^{a}(z)\bar{b}^{a}_{j\beta}(z)b_{j\alpha}^{a}(z)\Bigg{]},$
( 4.3)

where $z$ is the complex time variable along the SK contour $\mathcal{C}$
[Fig.4.1]. Here $Z$ is the equilibrium partition function which appears due to
the initial thermal density matrix Bera  _et al._ (2022). The initial density
matrix corresponds to the replica symmetric saddle point above $T_{SG}$, and
1RSB marginal SG saddle point below $T_{SG}$, as discussed in Sec.3.1. After
performing the disorder averaging, we can take the $n\rightarrow 0$ limit. The
resulting time-translation invariant dynamical correlations and responses are
identical to those obtained from the 1RSB saddle point solutions in the
marginal SG phase. Using this replicated generating function, the OTOC can now
be obtained as,

$\displaystyle F_{1}^{ab}(t_{1},t_{2})=$
$\displaystyle\frac{1}{N^{2}M^{2}}\sum_{ij,\alpha\beta}\langle yb^{\dagger
a}_{i\alpha 4}(t_{1})yb^{\dagger a}_{j\beta 3}(0)yb_{i\alpha
2}^{b}(t_{2})yb_{j\beta 1}^{b}(0)\rangle$ ( 4.4a) $\displaystyle
F_{2}^{ab}(t_{1},t_{2})=$
$\displaystyle\frac{1}{N^{2}M^{2}}\sum_{ij,\alpha\beta}\langle yb_{i\alpha
4}^{a}(t_{1})yb^{\dagger a}_{j\beta 3}(0)yb_{i\alpha 2}^{\dagger
b}(t_{2})yb^{b}_{j\beta 1}(0)\rangle$ ( 4.4b)

where the superscripts $1,2,3,4$ denote the branch on the SK contour and the
average $\langle\dots\rangle$ is with respect to replicated generating
function [Eq.(4.2)].

Figure 4.2:  Diagrammatic representation of the kernel equation in Eq.(4.6b)
for $\mathcal{O}(1/(NM))$ term, ($\mathcal{F}_{1,2}^{a}$) in the OTOC
$F_{1,2}^{a}(t_{1},t_{2})$. The solid horizontal lines denote dressed retarded
and advanced propagators $G_{R}$, $G_{A}$ and the vertical rung denotes the
Wightmann correlations $G_{lr}^{+}$ and $G_{lr}^{-}$. The dashed line
represents disorder averaging. We suppress replica indices of the vertices in
the ladder diagram to avoid cluttering.

The diagrams that contribute to the OTOC can be arranged in powers of
$1/(NM)$. $\mathcal{O}(1)$ diagrams are disconnected and do not contribute to
the growth of OTOC. The contribution, $\mathcal{F}_{\mu}^{aabb}$ at
$\mathcal{O}(1/(NM))$ form the ladder diagrams (Fig.4.2) grow exponentially,
i.e., for
$F_{\mu}^{aabb}(t_{1},t_{2})=F_{\mu}^{0aabb}(t_{1},t_{2})+(1/(NM))\mathcal{F}_{\mu}^{aabb}(t_{1},t_{2})+\mathcal{O}(1/(N^{2}M^{2}))$
($\mu=1,2$), $\mathcal{F}_{\mu}\sim e^{\lambda_{\mathrm{L}}t}$. The ladder
diagrams can be written in the form of a Bethe-Salpeter-like equation, similar
to SYK-type models Gu  _et al._ (2020); Maldacena and Stanford (2016);
Banerjee and Altman (2017),

$\displaystyle\mathcal{F}_{1}^{aabb}$
$\displaystyle(t_{1},t_{2})=\sum_{c}\int_{t_{3},t_{4}}\left[K_{11}^{aacc}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{1}^{ccbb}(t_{3},t_{4})\right.$
$\displaystyle\left.+K_{12}^{aacc}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{2}^{ccbb}(t_{3},t_{4})\right]$
( 4.5a) $\displaystyle\mathcal{F}_{2}^{aabb}$
$\displaystyle(t_{1},t_{2})=\sum_{c}\int_{t_{3},t_{4}}\left[K_{21}^{aacc}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{1}^{ccbb}(t_{3},t_{4})\right.$
$\displaystyle\left.+K_{22}^{aacc}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{2}^{ccbb}(t_{3},t_{4})\right].$
( 4.5b)

Here $\int_{t_{3},t_{4}}=\int dt_{3}dt_{4}$. In the chaotic growth regime,
$\lambda_{\mathrm{L}}^{-1}\lesssim t<\lambda_{\mathrm{L}}^{-1}\ln N$, the
propagators along the horizontal lines, from $t_{1}$ to $t_{3}$ and from
$t_{2}$ to $t_{4}$ in Fig.4.2 can be approximated by retarded and advanced
propagators Kitaev ; Maldacena and Stanford (2016); Gu and Kitaev (2019).
Since the retarded and advanced propagators are replica-diagonal in both PM
and SG phases, the kernel equations become,

$\displaystyle\mathcal{F}_{1}^{a}$ $\displaystyle(t_{1},t_{2})=\int
dt_{3}dt_{4}K_{11}^{a}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{1}^{a}(t_{3},t_{4})$
$\displaystyle+\int
dt_{3}dt_{4}K_{12}^{a}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{2}^{a}(t_{3},t_{4})$
( 4.6a) $\displaystyle\mathcal{F}_{2}^{a}$ $\displaystyle(t_{1},t_{2})=\int
dt_{3}dt_{4}K_{21}^{a}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{1}^{a}(t_{3},t_{4})$
$\displaystyle+\int
dt_{3}dt_{4}K_{22}^{a}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{2}^{a}(t_{3},t_{4}),$
( 4.6b)

where
$\mathcal{F}_{\mu}^{a}(t_{1},t_{2})\equiv\mathcal{F}_{\mu}^{aaaa}(t_{1},t_{2})$
and the kernels are given by

$\displaystyle K_{11}^{a}(t_{1},t_{2},t_{3},t_{4})=$ $\displaystyle
2\tilde{J}^{2}G_{A}(t_{31})G_{R}(t_{24})G_{lr}^{+}(t_{43})G_{lr}^{-}(t_{34})$
$\displaystyle K^{a}_{12}(t_{1},t_{2},t_{3},t_{4})=$
$\displaystyle\tilde{J}^{2}G_{A}(t_{31})G_{R}(t_{24})G_{lr}^{+}(t_{43})G_{lr}^{+}(t_{43})$
$\displaystyle K^{a}_{21}(t_{1},t_{2},t_{3},t_{4})=$
$\displaystyle\tilde{J}^{2}G_{R}(t_{13})G_{A}(t_{42})G_{lr}^{-}(t_{34})G_{lr}^{-}(t_{34})$
$\displaystyle K^{a}_{22}(t_{1},t_{2},t_{3},t_{4})=$ $\displaystyle
2\tilde{J}^{2}G_{R}(t_{13})G_{A}(t_{42})G_{lr}^{-}(t_{34})G_{lr}^{+}(t_{43}).$
( 4.7a)

Here all the Green’s functions are replica diagonal and $G_{lr}^{\pm}(t)$
refers to Wightmann correlations functions (Appendix A.2).

Using the exponential ansatz in the growth regime,
$\mathcal{F}_{\mu}^{a}(t_{1},t_{2})=e^{\lambda_{\mathrm{L}}(t_{1}+t_{2})/2}f_{\mu}^{a}(t_{1}-t_{2})$
and taking the Fourier transform of the kernel equations [Eq.(4.6)], leads to
an eigenvalue equations in the frequency space with eigenvalue,
$\lambda_{e}=1$. For the PM phase, these eigenvalue equations are

$\displaystyle\tilde{J}^{2}G_{A}\left(-\omega-\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)G_{R}\left(-\omega+\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)$
$\displaystyle\left[2\int
d\omega^{\prime}g_{1}(\omega-\omega^{\prime})\right.$
$\displaystyle\left.f_{1}^{a}(\omega^{\prime})+\int
d\omega^{\prime}g_{2}(-(\omega-\omega^{\prime}))f_{2}^{a}(\omega^{\prime})\right]$
$\displaystyle=f_{1}^{a}(\omega)$ ( 4.8a)
$\displaystyle\tilde{J}^{2}G_{R}\left(\omega+\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)G_{A}\left(\omega-\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)$
$\displaystyle\left[\int d\omega^{\prime}g_{2}(\omega-\omega^{\prime})\right.$
$\displaystyle\left.f_{1}^{a}(\omega^{\prime})+2\int
d\omega^{\prime}g_{1}(\omega-\omega^{\prime})f_{2}^{a}(\omega^{\prime})\right]$
$\displaystyle=f_{2}^{a}(\omega)$ ( 4.8b)

where

$\displaystyle g_{1}(\omega)$ $\displaystyle=\frac{1}{4\pi^{2}}\int
d\omega_{1}G_{lr}^{+}(\omega_{1})G_{lr}^{-}(\omega+\omega_{1})$
$\displaystyle=\frac{1}{4}\int
d\omega_{1}\frac{\rho(\omega_{1})\rho(\omega+\omega_{1})}{\mathrm{sinh}\frac{\beta\omega_{1}}{2}\mathrm{sinh}\frac{\beta(\omega+\omega_{1})}{2}}$
( 4.9a) $\displaystyle g_{2}(\omega)$ $\displaystyle=\frac{1}{4\pi^{2}}\int
d\omega_{1}G_{lr}^{+}(\omega_{1})G_{lr}^{+}(\omega-\omega_{1})$
$\displaystyle=\frac{1}{4}\int
d\omega_{1}\frac{\rho(\omega_{1})\rho(\omega-\omega_{1})}{\mathrm{sinh}\frac{\beta\omega_{1}}{2}\mathrm{sinh}\frac{\beta(\omega-\omega_{1})}{2}}.$
( 4.9b)

Here $\rho(\omega)=-(1/\pi)\mathrm{Im}G_{R}(\omega)$ is the bosonic spectral
function (Appendix A). We discretize the kernel equations [Eqs.(4.8)] in
frequency and diagonalize the resulting matrix kernel for different trial
values of $\lambda_{\mathrm{L}}$. Eventually, the Lyapunov exponent
$\lambda_{\mathrm{L}}$ is obtained from the one for which the kernel has at
least one eigenvalue $\lambda_{e}=1$.

Similarly, for the 1RSB SG phase, the kernel equations are

$\displaystyle\tilde{J}^{2}G_{A}\left(-\omega-\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)G_{R}\left(-\omega+\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)\bigg{[}2\Big{(}g^{2}f_{1}^{a}(\omega)+g\int\frac{d\omega^{\prime}}{2\pi}\tilde{g}_{lr}(\omega-\omega^{\prime})f_{1}^{a}(\omega^{\prime})+g\int\frac{d\omega^{\prime}}{2\pi}\tilde{g}_{lr}(\omega^{\prime}-\omega)f_{1}^{a}(\omega^{\prime})$
$\displaystyle+\int
d\omega^{\prime}\tilde{g}_{1}(\omega-\omega^{\prime})f_{1}^{a}(\omega^{\prime})\Big{)}+g^{2}f_{2}^{a}(\omega)+2g\int\frac{d\omega^{\prime}}{2\pi}\tilde{g}_{lr}(-\omega+\omega^{\prime})f_{2}^{a}(\omega^{\prime})+\int
d\omega^{\prime}\tilde{g}_{2}(-\omega+\omega^{\prime})f_{2}^{a}(\omega^{\prime})\bigg{]}=f_{1}^{a}(\omega)$
( 4.10a)
$\displaystyle\tilde{J}^{2}G_{A}\left(\omega-\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)G_{R}\left(\omega+\mathfrak{i}\frac{\lambda_{\mathrm{L}}}{2}\right)\bigg{[}2\Big{(}g^{2}f_{2}^{a}(\omega)+g\int\frac{d\omega^{\prime}}{2\pi}\tilde{g}_{lr}(\omega-\omega^{\prime})f_{2}^{a}(\omega^{\prime})+g\int\frac{d\omega^{\prime}}{2\pi}\tilde{g}_{lr}(\omega^{\prime}-\omega)f_{2}^{a}(\omega^{\prime})$
$\displaystyle+\int
d\omega^{\prime}\tilde{g}_{1}(\omega-\omega^{\prime})f_{2}^{a}(\omega^{\prime})\Big{)}+g^{2}f_{1}^{a}(\omega)+2g\int\frac{d\omega^{\prime}}{2\pi}\tilde{g}_{lr}(\omega-\omega^{\prime})f_{1}^{a}(\omega^{\prime})+\int
d\omega^{\prime}\tilde{g}_{2}(\omega-\omega^{\prime})f_{1}^{a}(\omega^{\prime})\bigg{]}=f_{2}^{a}(\omega)$
( 4.10b)

In the above equations, $\tilde{g}_{1}(\omega)$ and $\tilde{g}_{2}(\omega)$
are the same as in Eqs.(4.9),(4.9) in terms of the SG spectral function
$\rho(\omega)$. We obtain $\lambda_{\mathrm{L}}$ as a function of $S$ and $T$
by separately solving Eq.(4.8) for the PM phase, and Eq.(4.10)] for the SG
phase, with the retarded, advanced and Wightmann functions obtained from the
large-$N$ saddle point equations (Appendix A).

The crucial difference between the PM and SG phases is encoded in the ladder
kernerl [Eq.(4.7)] through the Wightmann correlators. In the SG phase,

$\displaystyle G_{lr}^{+}(\omega)=G_{lr}^{-}(\omega)=2\pi
g\delta(\omega)\epsilon_{ab}+\pi\rho(\omega)/(\mathrm{sinh}(\beta\omega/2))\delta_{ab}$
( 4.11)

In contrast, the delta function term is absent in the PM phase, where
$g=\sqrt{q_{\mathrm{EA}}}=0$.

Figure 4.3: Variation of Lyapunov exponent with $S$ in the SY model.
$\lambda_{\mathrm{L}}(S)$ for $T=0.1$, $T=0.8$ and $T=2.0$ (in units of
$\tilde{J}$), across the PM-SG phase transition marked by the vertical lines
at $S_{SG}\simeq 0.48,~{}1.29,~{}2.1$, respectively. $\lambda_{\mathrm{L}}$
has a maximum in the PM phase and changes non-monotonically.
$\lambda_{\mathrm{L}}\sim 1/S^{1.7}$ in the SG phase.

We compare the results for $\lambda_{\mathrm{L}}(S,T)$ obtained using the
above formalism for the zero-dimensional SY model with the
$\lambda_{\mathrm{L}}(\Gamma,T)$ for the effective zero-dimensional $p$-spin
glass model. The OTOC formalism for the latter is very similar to that of SY
model, but somewhat simpler. To extract $\lambda_{\mathrm{L}}(\Gamma,T)$, we
compute the following regularized OTOC Bera  _et al._ (2022) for the $p$-spin
glass model,

$\displaystyle F(t_{1},t_{2})=$
$\displaystyle\frac{1}{N^{2}}\sum_{i,j}\overline{\mathrm{Tr}\Big{[}ys_{i}(t_{1})ys_{j}(0)ys_{i}(t_{2})ys_{j}(0)\Big{]}},$
( 4.12)

We skip the details of the formalism here since those have been already
described in Ref.Bera  _et al._ , 2022. We briefly discuss a few additional
details for the non-local OTOC for $p$-spin glass chain model [Eq.(2.2)] in
Sec.5 along with those for the SY chain [Eq.(2.1)].

###  4.1 Lyapunov exponent across spin glass transition in the SY model

We obtain $\lambda_{\mathrm{L}}(S,T)$ in the PM and SG phases by separately
solving the corresponding large-$N$ saddle-point equations (Sec.3) for the two
phases. Specifically, as shown in Fig.3.1, we take cuts along several $S$ and
$T$, namely, from small $S$ (low-$T$) quantum critical region $a$. $S=0.5$
($d$. $T=0.1$), from intermediate-$S$ (intermediate-$T$) crossover region
between quantum critical and local moment behaviour $b$. $S=1$ ($e$. $T=0.8$),
and from large-$S$ (high-$T$) local moment region $c$. $S=2$ ($f$. $T=2$), as
a function of $T$ ($S$). We set $\tilde{J}=1$ as the unit of temperature
(energy). The results for $\lambda_{\mathrm{L}}(S,T)$ are qualitatively
similar for all the cuts, as we discuss below.

$S$ dependence of Lyapunov exponent.—We first show the Lyapunov exponent as a
function of $S$ in Fig.4.3 from PM to spin glass phase for three temperatures,
$T=0.1,0.8,2.0$. In all the cases, the Lyapunov exponent
$\lambda_{\mathrm{L}}$ has a non-monotonic dependence on $S$ within the PM
phase; $\lambda_{\mathrm{L}}$ initially increases with $S$ till a value
$S_{m}(T)$ where $\lambda_{\mathrm{L}}$ has a broad peak, and then it
decreases with further increase of $S$ approaching the critical value of $S$,
$S_{SG}(T)$ for the spin glass transition, which is obtained from the
breakpoint criterion $m(T,S_{SG})=1$. The non-monotonic
$\lambda_{\mathrm{L}}(S)$ in the PM phase can be heuristically understood from
the behaviour of the spectral function. As discussed in Appendix A.2, in our
numerical solution of the real-frequency saddle-point equation [Eq.(3.3)],
e.g., at $T=0.1,0.8$, for small $S\simeq 0.05-0.2$, the spectral function has
a local moment like form [see Figs.A.1(a), A.2(a)], with a peak at
$\omega\sim\omega_{0}\sim T\ln[(S+1)/S]>T$ and very little spectral weight
around $\omega=0$. The peak moves closer to zero energy with increasing $S$
and the spectral density increases around $\omega=0$ as the spectral function
resembles the spin-liquid solution [Eq.(A.14)] more. As a result,
$\lambda_{\mathrm{L}}$ is small at small $S$ and increases with $S$ till a
maximum value in the PM phase. However, beyond the maximum, closer to the SG
transition, there is a loss of spectral weight at intermediate energies or
_pseudogap_ Camjayi and Rozenberg (2003) like feature, in between a narrower
low-energy peak and a broad peak at high energies. This presumably leads to
decreases of $\lambda_{\mathrm{L}}$ approaching the SG transition. Thus,
$\lambda_{\mathrm{L}}$ exhibits a non-monotonic dependence on $S$ in the PM
region.

In the PM phase, especially in the quantum critical region $IV$ in Fig.3.1,
for $S\to 0$, we expect dominantly spin liquid behaviour and
$\lambda_{\mathrm{L}}(S\to 0,T)\to\lambda_{\mathrm{L}}^{\mathrm{SYK}}(T)$,
i.e., $\lambda_{\mathrm{L}}$ should asymptotically approach the $T$-dependent
Lyapunov exponent Banerjee and Altman (2017),
$\lambda_{\mathrm{L}}^{\mathrm{SYK}}(T)$, of the SYK spin liquid Georges  _et
al._ (2001) at $S=0$. Based on this expectation we fit our data with
$\lambda_{\mathrm{L}}(S,T)\simeq\lambda_{\mathrm{L}}(0,T)+aS^{\alpha_{\lambda}}$
in the PM phase for small values of $S$, before reaching the peak [Fig.4.3].
As shown in Fig.4.4(a) for $T=0.8,1.0$, we find reasonable fit with
$\alpha_{\lambda}\simeq 0.5-0.6$. However, due to this $\sqrt{S}$ dependence,
there is a rapid decrease of $\lambda_{\mathrm{L}}(S)$ for $S\lesssim 0.2$, as
seen in Fig.4.4(a). As a result, the extrapolated $\lambda_{\mathrm{L}}(0,T)$
for $S\to 0$ from the fit consistently turns out to be very low, and much
smaller than $\lambda_{\mathrm{L}}^{SYK}(T)$ Banerjee and Altman (2017), even
at low temperature $T\simeq 0.1$. Nevertheless, our real frequency saddle-
point solutions, e.g., the computed bosonic spectral function, for $S\leq
0.05$ at low $T$ are consistent with the conformal spin liquid solutions
Sachdev and Ye (1993); Georges  _et al._ (2001); Sachdev (2015) for
$|\omega|,T\ll J$ as well as the numerical spin liquid solutions obtained in
Ref. Camjayi and Rozenberg, 2003, as discussed in Appendix A.2. Thus, our
results indicate strong $S$ (and $T$) dependent corrections to the conformal
limit of the spectral function, and, consequently, $\lambda_{\mathrm{L}}=2\pi
T$, in the spin liquid region, even for small $S\lesssim 0.05$. We also note
that our numerical saddle-point solutions at small $S$ smoothly evolves to the
solutions at higher $S$, as evident from the smooth evolution of
$\lambda_{\mathrm{L}}(S)$ in Figs.4.3, 4.4(a) for all $T$.

Figure 4.4: $S$ dependence of $\lambda_{\mathrm{L}}$ in PM and SG phases of SY
model. $\lambda_{\mathrm{L}}$ is shown as a function of $S$. (a) In the PM
phase, before reaching the maximum [Fig.4.3], $\lambda_{\mathrm{L}}$ varies
as,
$\lambda_{\mathrm{L}}(S,T)\simeq\lambda_{\mathrm{L}}(0,T)+aS^{\alpha_{\lambda}}$,
shown as line fit to the data (circles) for $T=0.8,2.0$ with
$\alpha_{\lambda}\approx 0.5,0.6$, respectively. (b) $\lambda_{\mathrm{L}}$
decays with $S$ as a power law in the SG phase, shown as straight line fit to
the data (circles) on log-log scale. $\lambda_{\mathrm{L}}\sim 1/S^{1.7}$ for
$T=0.8,2.0$.

The Lyapunov exponent appears to have weak discontinuity/change of slope at
$S_{SG}$ for $T=0.1,~{}2.0$, which, however, could be a numerical artifact due
to problem in obtaining a converged solution of the saddle-point equations
very close to the transition, as well as, due to difficulty in matching the PM
and 1RSB SG solutions approaching from two different sides of the SG
transition. We find $\lambda_{\mathrm{L}}(S)$ to be continuous for $T=0.8$
across the SG transition. In the SG phase, $\lambda_{\mathrm{L}}(S)$ decreases
monotonically as $\sim 1/S^{1.7}$, as shown in Fig.4.4(b). The reduction of
$\lambda_{\mathrm{L}}$ with increasing $S$ or its increase with decreasing $S$
(increasing quantum fluctuations) can be understood qualitatively by drawing
analogy with the similar dependence of $\lambda_{\mathrm{L}}$ on $\Gamma$ in
the $p$-spin glass model. In the latter case, the kernel equations for
$\lambda_{\mathrm{L}}$ can be mapped to a solvable one-dimensional Schrödinger
equation Bera  _et al._ (2022). As a result, it can be shown by that
$\lambda_{\mathrm{L}}\sim T^{2}/q_{\mathrm{EA}}^{3/2}$. Thus
$\lambda_{\mathrm{L}}$ increases with increasing quantum fluctuations in the
SG phase due to the reduction of EA order parameter. Similar relation between
$\lambda_{\mathrm{L}}$ and $q_{\mathrm{EA}}$ can be expected for the SY model,
however, the kernel equations [Eq.(4.10)] are much more complicated in this
case, compared to the $p$-spin glass model Bera  _et al._ (2022), and we could
not map them into a solvable quantum mechanical problem.

Figure 4.5: Temperature dependence of Lyapunov exponent in SY model.
$\lambda_{\mathrm{L}}$ as function of $T$ for $S=0.5$, $S=1.0$ and $S=2.0$
across the SG to PM phase transitions, marked by the vertical lines at
$T_{SG}\simeq 0.12,~{}0.49,~{}1.9$, respectively. The inset shows zoomed in
view near the transitions for smaller $S$ values.

Temperature dependence of $\lambda_{\mathrm{L}}$.— We show the temperature
dependence of $\lambda_{\mathrm{L}}$ in Fig.4.5 for three values of $S$,
$S=0.5,~{}1.0,~{}2.0$. In all these case, in contrast to non-monotonic
dependence of $\lambda_{\mathrm{L}}$ on $S$, $\lambda_{\mathrm{L}}$
monotonically increases with temperature throughout SG and PM phases, with
possibly a weak slope change at the transition [Fig.4.5 (inset)]. The
dependence of $\lambda_{\mathrm{L}}$ on temperature, in SG and PM phases is
summarized below. Similar to the case of $p$-spin glass Bera  _et al._ (2022),
$\lambda_{\mathrm{L}}$ follows a power-law $T$ dependence in the SG phase for
SY model. For large $S$, $\lambda_{\mathrm{L}}\sim T^{\beta_{\lambda}}$, where
the exponent $\beta_{\lambda}\approx 1.5$ weakly depends on $S$, as shown in
Fig.4.6. For small $S$, $\lambda_{\mathrm{L}}\sim T$, i.e., has a linear $T$
dependence [Fig.4.7(a)], albeit in a narrow range of temperature accessible
for $T\lesssim T_{SG}$ for small $S$.

In the PM phase for $T>T_{SG}$, for all values of $S$, $\lambda_{\mathrm{L}}$
increases near the transition logarithmically in $T$ as, $a\ln{T}+c$, with
constants $a$ and $c$ that depend on $S$ [Fig.4.7(b)].
$\lambda_{\mathrm{L}}(T)$ eventually tends to saturate at high temperature, as
shown in Fig.4.7. The logarithmic dependence of $\lambda_{\mathrm{L}}(T)$,
albeit over a small temperature window near $T_{SG}$ in the PM phase, even for
relatively large $S=2$, is somewhat surprising. For classical spin models
Bilitewski  _et al._ (2018); Ruidas and Banerjee (2021),
$\lambda_{\mathrm{L}}\sim\sqrt{T}$ is found at high temperatures, e.g., above
thermal phase transitions. The $\sqrt{T}$ behaviour is expected for classical
models from fairly generally considerations Kurchan (2016). In our case, the
initial increase of $\lambda_{\mathrm{L}}(T)$ for $T\gtrsim T_{SG}$ cannot be
fitted with a power law like $\sqrt{T}$. Nevertheless, qualitatively, the
increase of $\lambda_{\mathrm{L}}$ is expected with increasing $T$ due to
gradual unfreezing of spins going away from $T_{SG}$. As discussed in Appendix
A.2, near $T_{SG}$ in the PM phase, the spectral function has a sharp peak
close to $\omega=0$ separated from high-frequency spectral weight with a
pseudogap-like dip at intermediate frequencies. The sharp peak can be
understood as a relic of the delta function peak $\sim\delta(\omega)$ that
appears in correlation functions (see Appendix A.2) in the SG phase, e.g., in
the Wightmann correlation function [Eq.(4.11)]. The delta function peak
indicates static order and spin freezing. With increasing $T$, the low- and
high-frequency peaks move closer filling up the pseudogap at intermediate
frequencies due to unfreezing of spins. This spectral weight transfer
correlates with increase of $\lambda_{\mathrm{L}}$ with $T$. Like in the SYK
model Banerjee and Altman (2017), $\lambda_{\mathrm{L}}$ is expected to
saturate eventually at high temperature since the model has a bounded
spectrum. This is unlike the $p$-spin glass model which has a unbounded
spectrum due to unbounded kinetic energy. The latter renders interaction
effect negligible Bera  _et al._ (2022) at high temperature in the $p$-spin
glass model leading to decrease of $\lambda_{\mathrm{L}}\sim 1/T^{2}$ for
$T\gg\tilde{J}$.

Figure 4.6:  Power-law temperature dependence of $\lambda_{\mathrm{L}}$ in the
SG phase for large $S$ in SY model. The straight lines show the power-law fit,
$\lambda_{\mathrm{L}}\sim T^{\beta_{\lambda}}$ to the data (circles), on log-
log plot with $\beta_{\lambda}=1.6,1.5$ for $S=2.0,3.0$, respectively. Figure
4.7: Temperature dependence of $\lambda_{\mathrm{L}}$ in PM and SG phases of
SY model. (a) $\lambda_{\mathrm{L}}$ increases linearly with temperature as
$\lambda_{\mathrm{L}}\sim T$ in the SG phase for small $S$, shown for
$S=0.5,1.0$ as lines fitted to the data (circles) on log-log scale. (b) In the
PM phase, $\lambda_{\mathrm{L}}$ shows logarithmic dependence on temperature,
$a\ln{T}+c$ (line) initially before it saturates to a constant value for large
$T$ (Fig.4.5), as shown for $S=0.5,1.0,2.0$ in semi-log plot
($\lambda_{\mathrm{L}}$ vs. $\ln{T}$). Figure 4.8: Variation of Lyapunov
exponent with quantum parameter $\Gamma$ in the $p$-spin glass model.
$\lambda_{\mathrm{L}}$ as function of $\Gamma$ for $T=0.2$ and $T=0.4$ across
the PM-SG phase transition marked by the vertical lines at $\Gamma_{SG}\simeq
2.65,~{}1.35$ respectively. $\lambda_{\mathrm{L}}$ changes non-monotonically
and has a maximum above the the glass transition in the PM phase Bera  _et
al._ (2022). Figure 4.9: Temperature dependence of $\lambda_{\mathrm{L}}$ in
the $p$-spin glass model. $\lambda_{\mathrm{L}}$ as function of $T$ for
$\Gamma=0.5$ and $\Gamma=2.0$ across the PM-SG phase transition marked by the
vertical lines at $T_{SG}\sim 0.3,0.51$ respectively. $\lambda_{\mathrm{L}}$
changes non-monotonically and has a maximum above the glass transition in the
PM phase Bera  _et al._ (2022)

The non-monotonic behaviour of $\lambda_{\mathrm{L}}(S)$ [Fig.4.5] in the PM
phase of the SY model is similar to the non-monotonic behaviour seen Bera  _et
al._ (2022) for $\lambda_{\mathrm{L}}(\Gamma)$ or
$\lambda_{\mathrm{L}}(\hbar)$, i.e., as a function of quantum parameter, in
the $p$-spin glass model [Eq.(2.2)] as shown in Fig.4.8. However, in the
latter model, in contrast to the $\lambda_{\mathrm{L}}(T)$ [Fig.4.3] in the SY
model, $\lambda_{\mathrm{L}}$ also shows non-monotonic dependence on $T$ in
the PM phase, as in Fig.4.9. The broad peak related to the non-monotonic
$\lambda_{\mathrm{L}}(\Gamma,T)$ in the PM phase of the $p$-spin glass model
correlates with the onset of complex glassy relaxation Bera  _et al._ (2022),
where the dynamics starts getting dominated by saddle points of the glassy
energy landscape Correale  _et al._ (2023). In the classical limit ($\Gamma\to
0$) of the $p$-spin glass model, the origin of the broad peak was shown to be
related with the interplay of the rapid increase of relaxation time
approaching the glass transition, and weakening of interaction effects due to
unbounded kinetic energy growth with increasing temperature. The interplay
leads to a crossover from strong to weak chaos with increasing $T$, manifested
as a broad peak in $\lambda_{\mathrm{L}}$ coinciding with the onset of two-
step glassy relaxation regime Bera  _et al._ (2022). Due to bounded spectrum
of the SY model, such crossover from strong to weak chaos with temperature is
absent in the SY model, and the $\lambda_{\mathrm{L}}(T)$ monotonically
increases with $T$ in the PM phase till it saturates at high temperature
[Fig.4.3]. However, as already discussed, a non-monotonic
$\lambda_{\mathrm{L}}(S)$ is still observed in the PM phase due to change of
spectral properties as a function of $S$.

Figure 4.10: Onset of glassy relaxation with $S$ in SY model. Correlation
function $C(t)$ in the PM phase near the PM-SG phase transition for a fixed
value $T=0.8$ and different values of $S$. As we approach the transition from
the PM phase, near the transition at $S_{SG}\simeq 1.3$, the two step glassy
relaxation is more pronounced. Figure 4.11: Onset of glassy relaxation with
temperature in SY model. Correlation function $C(t)$ in the PM phase near the
PM-SG phase transition for a fixed value of $S=1.0$ and different values of
$T$. The two step glassy relaxation is more pronounced as we approach
$T_{SG}=0.5$ from the PM phase.

The two-step glassy relaxation, seen for the $p$-spin glass model Bera  _et
al._ (2022), can also be observed in the PM phase of the SY model above the
glass transition. To this end, we calculate the correlation function using the
Fluctuation dissipation theorem,
$C(t)=-(\mathfrak{i}/M)\sum_{\alpha}\langle[b_{i\alpha
x}(t),b^{\dagger}_{i\alpha
x}(0)]\rangle=2\int(d\omega/2\pi)e^{-\mathfrak{i}\omega
t}\mathrm{coth}(\beta\omega/2)\mathrm{Im}G_{R}(\omega)$. We show the
correlation function for several $S$ in Fig .4.10 at $T=0.8$, and for several
$T$ in Fig .4.11 at $S=1$, above the SG transition. The gradual onset of the
two-step relaxation is evident. This suggests that the onset of two-step
glassy relaxation is not a sufficient condition to realize a peak in the
Lyapunov exponent.

##  5 Butterfly velocity in quantum spin glass chains

In this section we study the spatial spreading of chaos, first in the SY model
and then in the $p$-spin glass model.

Figure 5.1:  Diagrammatic representation of the Kernel equation in Eq.5.2 for
$\mathcal{O}(1/(NM))$ term, $\mathcal{F}_{\mu,x}^{a}(t_{1},t_{2})$ in the OTOC
$F_{\mu,x}^{a}(t_{1},t_{2})$ ($\mu=1,2$). The solid horizontal lines denote
dressed retarded and advanced propagators $G_{R}$, $G_{A}$ and the vertical
rung denotes the Wightmann correlations $G_{lr}^{+}$ and $G_{lr}^{-}$. The
dashed and dotted lines denote disorder averaging over $J$ and $J^{\prime}$
couplings, respectively in Eq.(2.1). Here $x^{\prime}=x\pm 1$. We suppress
replica indices for the vertices to avoid cluttering.

###  5.1 Butterfly velocity across spin glass transition in the SY chain

In the previous section, we calculated $\lambda_{\mathrm{L}}$ which determines
the temporal growth of the four-point OTOC. Due to the spatial locality of our
generalized 1D SY chain model [Eq.(2.1)], we can study the spatio-temporal
dynamics of growth and spread of chaos. The quantity which measures the
propensity of spatial spreading of chaos is the butterfly velocity, $v_{B}$.
The operators that form large commutators with an initial operator lie inside
a region of space-time that forms a _light cone_ Gu  _et al._ (2017); Khemani
_et al._ (2018), signifying a ballistic spread of chaos. The speed or velocity
that defines the boundary of the light cone is defined as $v_{B}$. Like the
Lyapunov exponent $\lambda_{\mathrm{L}}$, the butterfly speed $v_{B}$ provides
a state-dependent measure of operator growth and spreading, and can be related
Roberts and Swingle (2016) to the well-known Lieb-Robinson velocity bound Lieb
and Robinson (1972), that gives state-independent measure in terms of operator
norm. To this end, we calculate the following disorder-averaged regularized
spatio-temporal OTOCs for the SY model,

$\displaystyle F_{1,x}(t_{1},t_{2})$
$\displaystyle=\frac{1}{N^{2}M^{2}}\sum_{ij\alpha\beta}\overline{\mathrm{Tr}[yb^{\dagger}_{i\alpha,x}(t_{1})yb^{\dagger}_{j\beta,0}(0)yb_{i\alpha,x}(t_{2})yb_{j\beta,0}(0)]}$
( 5.1a) $\displaystyle F_{2,x}(t_{1},t_{2})$
$\displaystyle=\frac{1}{N^{2}M^{2}}\sum_{ij\alpha\beta}\overline{\mathrm{Tr}[yb_{i\alpha,x}(t_{1})yb^{\dagger}_{j\beta,0}(0)yb_{i\alpha,x}^{\dagger}(t_{2})yb_{j\beta,0}(0)]}$
( 5.1b)

As discussed in Sec.4 for the case of the onsite OTOCs, which is equivalent to
OTOCs in zero-dimensional SY model, the above non-local OTOCs are computed via
a replicated Keldysh field theory, where the OTOCs become replica diagonal up
to $\mathcal{O}(1/(NM))$ diagrams. The latter are arranged as ladder series
that can be written in the form of Bethe-Salpeter-like equations, as shown in
Fig:5.1,

$\displaystyle\mathcal{F}_{1,x}^{a}(t_{1},t_{2})=\int dx^{\prime}\int
dt_{3}dt_{4}\bigg{[}K_{11,xx^{\prime}}^{a}(t_{1},t_{2},t_{3},t_{4})$
$\displaystyle~{}\mathcal{F}_{1,x^{\prime}}^{a}(t_{3},t_{4})+K_{12,xx^{\prime}}^{a}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{2,x^{\prime}}^{a}(t_{3},t_{4})\bigg{]},$
( 5.2a) $\displaystyle\mathcal{F}_{2,x}^{a}(t_{1},t_{2})=\int dx^{\prime}\int
dt_{3}dt_{4}\bigg{[}K_{21,xx^{\prime}}^{a}(t_{1},t_{2},t_{3},t_{4})$
$\displaystyle~{}\mathcal{F}_{1,x^{\prime}}^{a}(t_{3},t_{4})+K_{22,xx^{\prime}}^{a}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}_{2,x^{\prime}}^{a}(t_{3},t_{4})\bigg{]}.$
( 5.2b)

In the above, the kernels, which are diagrammatically shown in Fig.5.1, are
given in Eqs.(5.3). Although we have not indicated explicitly, all the Green’s
functions appearing in the kernel elements are replica diagonal. Since the
saddle-point Green’s function (Sec.3.1) of the extended 1D SY chain model are
local and independent of $x$, the kernel elements are translation invariant,
as in Eqs.(5.3). Thus, after Fourier transforming Eqs.(5.2) to momentum space,
we obtain equations for the Fourier transform of
$\mathcal{F}^{a}_{\mu,x}(t_{1},t_{2})$ ($\mu=1,2$),
$\mathcal{F}^{a}_{\mu,p}(t_{1},t_{2})$. The kernel elements for a given
momentum $p$ are obtained by replacing $\tilde{J}^{2}$ in the local kernel
elements of Eqs.(4.7) by momentum dependent prefactors, as shown in Eqs.(5.4),
where $\alpha=J^{\prime 2}/\tilde{J}^{2}$. In this section, as discussed
below, we calculate butterfly velocity $v_{B}$ using two methods –(1) a fully
numerical method, which uses a straightforward generalization of the
exponential growth ansatz for the local OTOC (Sec.4) for each momentum mode;
we refer to this as _real $p$ method_, and (2) a semi-analytic method, called
the _single mode ansatz_ , following ref.Gu  _et al._ , 2017. Here we refer to
this as _imaginary $p$ method_.

$\displaystyle
K_{11,xx^{\prime}}^{a}(t_{1},t_{2};t_{3},t_{4})=G_{A}(t_{31})G_{R}(t_{24})G_{lr}^{+}(t_{43})G_{lr}^{-}(t_{34})\Big{(}2J^{2}\delta_{x,x^{\prime}}+J^{\prime
2}(\delta_{x,x^{\prime}}+\delta_{x\pm 1,x^{\prime}})\Big{)}$ ( 5.3a)
$\displaystyle
K_{12,xx^{\prime}}^{a}(t_{1},t_{2};t_{3},t_{4})=G_{A}(t_{31})G_{R}(t_{24})G_{lr}^{+}(t_{43})G_{lr}^{+}(t_{43})\Big{(}J^{2}\delta_{x,x^{\prime}}+J^{\prime
2}(\delta_{x\pm 1,x^{\prime}})\Big{)}$ ( 5.3b) $\displaystyle
K_{21,xx^{\prime}}^{a}(t_{1},t_{2};t_{3},t_{4})=G_{R}(t_{13})G_{A}(t_{42})G_{lr}^{-}(t_{34})G_{lr}^{-}(t_{34})\Big{(}J^{2}\delta_{x,x^{\prime}}+J^{\prime
2}(\delta_{x\pm 1,x^{\prime}})\Big{)}$ ( 5.3c) $\displaystyle
K_{22,xx^{\prime}}^{a}(t_{1},t_{2};t_{3},t_{4})=G_{R}(t_{13})G_{A}(t_{42})G_{lr}^{-}(t_{34})G_{lr}^{+}(t_{43})\Big{(}2J^{2}\delta_{x,x^{\prime}}+J^{\prime
2}(\delta_{x,x^{\prime}}+\delta_{x\pm 1,x^{\prime}})\Big{)}$ ( 5.3d)
$\displaystyle K_{11,p}^{a}(t_{1},t_{2};t_{3},t_{4})$
$\displaystyle=\tilde{J}^{2}\Big{(}2-\alpha(1-2\text{ }\mathrm{cos}\text{
}p)\Big{)}G_{A}(t_{31})G_{R}(t_{24})G_{lr}^{+}(t_{43})G_{lr}^{-}(t_{34})$ (
5.4a) $\displaystyle K_{12,p}^{a}(t_{1},t_{2};t_{3},t_{4})$
$\displaystyle=\tilde{J}^{2}\Big{(}1-\alpha(1-2\text{ }\mathrm{cos}\text{
}p)\Big{)}G_{A}(t_{31})G_{R}(t_{24})G_{lr}^{+}(t_{43})G_{lr}^{+}(t_{43})$ (
5.4b) $\displaystyle K_{21,p}^{a}(t_{1},t_{2};t_{3},t_{4})$
$\displaystyle=\tilde{J}^{2}\Big{(}1-\alpha(1-2\text{ }\mathrm{cos}\text{
}p)\Big{)}G_{R}(t_{13})G_{A}(t_{42})G_{lr}^{-}(t_{34})G_{lr}^{-}(t_{34})$ (
5.4c) $\displaystyle K_{22,p}^{a}(t_{1},t_{2};t_{3},t_{4})$
$\displaystyle=\tilde{J}^{2}\Big{(}2-\alpha(1-2\text{ }\mathrm{cos}\text{
}p)\Big{)}G_{R}(t_{13})G_{A}(t_{42})G_{lr}^{-}(t_{34})G_{lr}^{+}(t_{43})$ (
5.4d)

Figure 5.2: The growing part of OTOC, $\mathcal{F}_{1,x}(t)$, is plotted on
the $x-t$ plane in logarithmic scale, as indicated in the color bar at $T=0.8$
and $S=1.0$ using the real $p$ method.

Real $p$ method.— As in Sec.4, we take an exponential growth ansatz,
$\mathcal{F}_{\mu,p}(t_{1},t_{2})=e^{\lambda_{\mathrm{L}}(p)(t_{1}+t_{2})/2}f_{\mu,p}^{a}(t_{1}-t_{2})$,
and extract the momentum-dependent Lyapunov exponent $\lambda_{\mathrm{L}}(p)$
by demanding that the matrix kernel has at least one eigenvalue equal to 1
(typically degenerate). We compute the eigenvector corresponding to the
eigenvalue 1 to obtain $f_{\mu,p}^{a}(\omega)$. To extract $v_{B}$, we
calculate
$\mathcal{F}^{a}_{\mu,p}(t,t)=e^{\lambda_{\mathrm{L}}(p)t}f^{a}_{\mu,p}(0)$,
where $f_{\mu,p}(0)=\int(d\omega/2\pi)f_{\mu,p}(\omega)$. We Fourier transform
to the position space,
$\mathcal{F}_{\mu,x}(t)\equiv\mathcal{F}^{a}_{\mu,x}(t,t)=(1/L)\sum_{p}e^{ixp}\mathcal{F}_{\mu,p}(t,t)$
where $L$ is the number of sites or the number of discrete $p$ values taken
within the Brillouin zone $[-\pi,\pi]$ for unit lattice spacing. We show in
Fig.5.2 that $\mathcal{F}_{1,x}(t)$ exhibits quite distinctive light cone. The
magnitude of $\mathcal{F}_{\mu,p}(t,t)$ cannot be fixed, i.e., the eigenvector
can be multiplied by an arbitrary factor, for the kernel equations. Thus, we
compute the light-cones using the locus of $(x,t)$ points corresponding to
$\mathcal{F}_{1,x}(t)=\mathcal{F}_{th}$, for various threshold values
$\mathcal{F}_{th}$, for which $t\propto x$. We find that this ballistic
feature is violated for too small or too large values of $\mathcal{F}_{th}$.
We extract $v_{B}(\mathcal{F}_{th})$ from the resulting light cones and find
that $v_{B}(\mathcal{F}_{th})$ converges approximately to the same value
$v_{B}$ over a wide range of $\mathcal{F}_{th}$, as shown in Fig.5.3. We
discuss the results for $v_{B}(S,T)$ obtained from the real $p$ method below
along with that from the imaginary $p$ method.

Figure 5.3:  Light cones at $T=0.8$ and $S=1.0$ for different values of
$\mathcal{F}_{th}$ obtained using the real $p$ method. The same butterfly
speed $v_{B}$ is obtained for a wide range of $\mathcal{F}_{th}$.

Imaginary $p$ method.— To calculate butterfly velocity through this method, we
use the single-mode ansatz for OTOC and the results obtained by Gu and Kitaev
Gu and Kitaev (2019). Using the single mode ansatz Kitaev and Suh (2018) for
early time OTOC, they proved a ladder identity which implies the existence of
a pole and a saddle at two imaginary values of momentum that contribute to the
growth of OTOC. We refer to ref.Gu and Kitaev, 2019 for the proof of the
ladder identity and the details of the method. Using this ansatz we can write
the solution of Eqs.(5.2) for the 1D SY chain model [Eq.(2.1)],

$\displaystyle\mathcal{F}_{\mu,p}^{a}(t_{1},t_{2})\approx\frac{e^{\lambda_{\mathrm{L}}(p)(t_{1}+t_{2})/2}}{C(p)}\Upsilon_{\mu,p}^{R}(t_{12})\Upsilon_{\mu,p}^{A}(0)$
( 5.5)

where $\Upsilon_{\mu,p}^{R}(t)$ and $\Upsilon_{\mu,p}^{A}(t)$ are the retarded
and advanced vertex functions Gu and Kitaev (2019) which modify the OTOC by an
overall magnitude. The dominant $p$ dependence of the above comes from a pole
$\lambda_{\mathrm{L}}=2\pi T$ ($k_{\mathrm{B}}=1$, $\hbar=1$) in $C(p)$, since
$C(p)\propto\cos{[\lambda_{\mathrm{L}}(p)/4T]}$ Gu and Kitaev (2019). This
pole leads to the maximal chaos for non-Fermi liquids in SYK and related
models Gu  _et al._ (2017); Gu and Kitaev (2019); Guo  _et al._ (2019). The
asymptotic form of Eq.(5.5) ignores the initial correlations and any non-
linear effects Gu and Kitaev (2019); Kitaev and Suh (2018) and is valid only
at the butterfly or chaos front, i.e.,
$1/N\ll\mathcal{F}_{\mu,p}(t_{1},t_{2})/N\ll 1$. Fourier transforming Eq.(5.5)
back to real space, we obtain

$\displaystyle\mathcal{F}_{\mu,x}(t)\equiv\mathcal{F}^{a}_{\mu,x}(t,t)\approx\int_{-\infty}^{\infty}\frac{dp}{2\pi}\frac{e^{\lambda_{\mathrm{L}}(p)t+ipx}}{\cos(\lambda_{\mathrm{L}}(p)/4T)}f(p),$
( 5.6)

where $f(p)$, which depends on the kernel [Eqs.(5.4)] and vertex functions in
Eq.(5.5), is assumed to be analytical function of $p$ and non vanishing in the
momentum regime of interest Gu and Kitaev (2019).

As in the real $p$ method discussed above, we choose a fixed, albeit
arbitrary, value of $\mathcal{F}_{\mu,x}(t,t)\approx 1$, which is at the lower
end of the butterfly front Gu and Kitaev (2019), to obtain the light cone from
the resulting locus of $(x,t)$ points and extract $v_{B}$. For large $x$ and
$t$, the integral in Eq.(5.6) can be estimated by saddle-point method. The
saddle-point $p_{s}$ for the integrand is obtained from
$\lambda^{\prime}_{\mathrm{L}}(p_{s})t+\mathfrak{i}x=0$, where
$\lambda_{\mathrm{L}}^{\prime}(p)$ is the derivative. This together with the
condition $\mathcal{F}_{\mu,x}(t,t)\approx 1$, implying
$\lambda_{\mathrm{L}}(p)t+\mathfrak{i}px=0$, leads to the following equation
for $p_{s}$

$\displaystyle\lambda^{\prime}_{\mathrm{L}}(p_{s})=\frac{\lambda_{\mathrm{L}}(p_{s})}{p_{s}},$
( 5.7)

and a corresponding velocity

$\displaystyle
v_{s}=\left(\frac{x}{t}\right)_{s}=\mathfrak{i}\lambda^{\prime}_{\mathrm{L}}(p_{s}).$
( 5.8)

The saddle point leads to a contribution
$\mathcal{F}_{\mu,x}^{s}(t)\sim\exp{[\lambda_{\mathrm{L}}(p_{s})t+\mathfrak{i}p_{s}x]}$.
Similarly, the pole $p_{1}$ of the integrand in Eq.(5.6), such that
$\lambda_{\mathrm{L}}(p_{1})=2\pi T$, gives a contribution
$\mathcal{F}_{\mu,x}^{p}(t)\sim\exp{[2\pi Tt+\mathfrak{i}p_{1}x]}$ with an
associated velocity

$\displaystyle v_{m}=\mathfrak{i}\frac{2\pi T}{p_{1}}.$ ( 5.9)

To give rise to a well defined light cone, $\mathcal{F}_{\mu,x}^{s}(t)$
($\mathcal{F}_{\mu,x}^{p}(t)$) should be non-oscillatory in space and time,
implying that $p_{s}$ ($p_{1}$) is purely imaginary, i.e.,
$p_{s}=\mathfrak{i}|p_{s}|$ ($p_{1}=\mathfrak{i}|p_{1}|$), and
$\lambda_{\mathrm{L}}(\mathfrak{i}|p_{s}|)$ is purely real. In Appendix C, we
show that this is indeed the case by solving for $\lambda_{\mathrm{L}}(p)$ for
purely imaginary $p=\mathfrak{i}|p|$ from Eqs.(5.2) using the exponential
growth ansatz, as discussed earlier for the real $p$ method.
$\lambda_{\mathrm{L}}(\mathfrak{i}|p|)$ is found to be a convex function of
$|p|$ as expected Gu and Kitaev (2019). The light cone and the associated
$v_{B}$ are determined by either the saddle-point or the pole contribution to
$\mathcal{F}_{\mu,x}(t)$ depending on whether $|p_{s}|<|p_{1}|$ or
$|p_{s}|>|p_{1}|$, respectively. The conditions also translate to
$v_{s}(\equiv v(p_{s}))<v_{1}(\equiv v(p_{1}))$ or $v_{s}>v_{1}$,
respectively, where the velocity
$v(p)=\mathfrak{i}\lambda^{\prime}_{\mathrm{L}}(p)$, since
$\lambda_{\mathrm{L}}(\mathfrak{i}|p|)$ is a convex (up) function of $p$
(Appendix C).

_For the condition $|p_{s}|<|p_{1}|$_, the deformed contour, that goes through
the saddle point $\mathfrak{i}|p_{s}|$ on the imaginary axis in the evaluation
of the integral in Eq.(5.6), does not enclose the pole and hence the latter
does not contribute. As discussed in refs.Gu and Kitaev, 2019; Guo  _et al._ ,
2019, this condition is realized when $\lambda_{\mathrm{L}}(0)\ll 2\pi T$,
i.e., when the chaos is far away from maximal. In this case
$\mathcal{F}_{\mu,x}(t)\approx\mathcal{F}^{s}_{\mu,x}(t)$ and $v_{B}=v_{s}$.
By making a small $p$ approximation,
$\lambda_{\mathrm{L}}(p)\simeq\lambda_{\mathrm{L}}(0)-\tilde{D}p^{2}$ (with
$\tilde{D}>0$), we obtain
$v_{B}=v_{s}=[4\tilde{D}\lambda_{\mathrm{L}}(0)]^{1/2}$ and

$\displaystyle\mathcal{F}_{\mu,x}(t)\sim
e^{\lambda_{\mathrm{L}}(0)t(1-x^{2}/v_{B}^{2}t^{2})}$ ( 5.10)

Here the quantity $\tilde{D}=v_{B}^{2}/4\lambda_{\mathrm{L}}(0)$ has the
dimension of diffusion constant. For certain holographic theories Blake (2016)
and strongly interacting diffusive metals built from SYK-type quantum dots Gu
_et al._ (2017); Guo  _et al._ (2019), $4\tilde{D}$ has been found to be
exactly equal to energy diffusion constant at low temperature. $\tilde{D}$ has
also been found to be same as spin diffusion constant at high temperature in a
classical spin liquid Bilitewski  _et al._ (2018), and even for classical XY
model Ruidas and Banerjee (2021). Assuming that such relation between
$\tilde{D}$ and actual diffusion constant hold, at least, over a limited range
of $(S,T)$ in the phased diagram of SY chain, we use $\tilde{D}$ as a proxy
for diffusion constant and plot it as a function of $S$ and $T$ below.

_For the other condition $|p_{s}|>|p_{1}|$_, realized for
$\lambda_{\mathrm{L}}(0)\approx 2\pi T$, the pole contributes and dominates
over the saddle-point contribution, such that

$\displaystyle\mathcal{F}_{\mu,x}(t)\approx\mathcal{F}^{p}_{\mu,x}(t)\sim
e^{2\pi Tt(1-|x|/v_{B}t)}$ ( 5.11)

and $v_{B}=v_{m}$, as can be deduced from $\mathcal{F}_{\mu,x}(t)\approx 1$.
The above corresponds to maximal chaos, realized for SYK-type non-Fermi
liquids Gu  _et al._ (2017); Gu and Kitaev (2019); Guo  _et al._ (2019).

For most part of the phase diagram [Fig.3.1] we have studied, we see that
saddle point, $|p_{s}|$ remains well below the pole, $|p_{1}|$ and so the OTOC
receives dominant contribution from the saddle point. See Appendix C for more
details. Therefore the chaos front travels with a speed
$v_{B}=\mathfrak{i}\lambda_{\mathrm{L}}^{\prime}(\mathfrak{i}|p_{s}|)$ and the
Lyapunov exponent $\lambda_{\mathrm{L}}$ is non-maximal, as was verified
through direct calculation of $\lambda_{\mathrm{L}}(S,T)$ in Sec.4. The SY
chain can be analytically shown to be maximally chaotic in the limit $S\to
0,~{}T\to 0$ following the same methods Kitaev ; Maldacena and Stanford
(2016); Banerjee and Altman (2017) used for calculating
$\lambda_{\mathrm{L}}(T)$ in the SYK model.

Figure 5.4:  Variation of butterfly velocity with $S$ in the SY model. $v_{B}$
across the SG phase transition ($S_{SG}$, vertical dashed lines) as a function
of $S$ at $T=0.1,0.8$ (in units of $\tilde{J}=1$) and $\alpha=0.25$,
calculated using both real and imaginary $p$ methods, as indicated in the
legends. For the real $p$ method, we have used $\mathcal{F}_{th}=0.1$ for
$T=0.1$ and $\mathcal{F}_{th}=0.5$ for $T=0.8$, with ranges (error bars)
$\mathcal{F}_{th}=0.5-2.0$ and $\mathcal{F}_{th}=0.1-1.0$, respectively, for a
lattice of $L=50$ sites. Figure 5.5:  $S$ dependence of $v_{B}$ in the PM and
SG phases of SY model. (a) The The data (circles) for $v_{B}(S,T)$ in the PM
phase, before it reaches the maximum [Fig.5.4], is fitted with
$v_{B}(0,T)+aS^{\alpha_{v}}$ (line) for $T=0.8$ with $\alpha_{v}\approx 0.5$.
(b) Exponential dependence of $v_{B}$ on $S$ in the SG phase is shown as
straight line on a semi-log plot ($\ln{v_{B}}$ vs. $S$) for $T=0.8$.

We calculate $v_{B}$ as a function of $T$ and $S$, across the SG-PM phase
transition, for $S=0.5,1.0$ and $T=0.1,0.8$ along the $a,~{}b$ and $d,~{}e$
cuts, respectively, shown in Fig.3.1, using both real and imaginary $p$
methods. As discussed below, we find that the results from the two methods
match closely, validating the particular form of single-mode ansatz in
Eq.(5.5) even in the replica-symmetry broken spin glass phase.

Figure 5.6: $\tilde{D}\sim v_{B}^{2}/\lambda_{\mathrm{L}}$ across the phase
transition ($S_{SG}$, vertical dashed lines) as a function of $S$ calculated
for $T=0.1$ and $T=0.8$.

$S$ dependence of butterfly velocity.— Fig.5.4 shows $v_{B}$ as a function of
$S$ for $T=0.1$ and $T=0.8$, going from spin liquid and spin liquid-local
moment crossover region to the SG phase. As in the case of
$\lambda_{\mathrm{L}}(S)$ [Fig.4.3], $v_{B}(S)$ exhibits a broad maximum in
the PM phase. However, unlike $\lambda_{\mathrm{L}}(S)$ it has a dip at the
phase transition $S_{SG}$ and then it increases monotonically with $S$, as the
system becomes more classical. Such dip in $v_{B}$ at the phase transition has
been also seen for classical spin models Ruidas and Banerjee (2021);
Bilitewski  _et al._ (2021), implying that the speed of information scrambling
becomes slower at the transition in all these cases.

Figure 5.7: Temperature dependence of butterfly velocity in SY model. $v_{B}$
across the phase transition ($T_{SG}$, vertical dashed lines) as a function of
temperature at $S=0.5,1.0$, $\alpha=0.25$, using both real and imaginary $p$
methods (see main text). The parameters and the conventions are same as those
described in the caption of Fig.5.4 .

We find $v_{B}$ to have power-law dependence,
$v_{B}(S,T)-v_{B}(0,T)\sim\sqrt{S}$ for small $S$, similar to
$\lambda_{\mathrm{L}}$ [Fig.4.4(a)] with a very small extrapolated
$v_{B}(0,T)\approx 0$ in the limit $S\to 0$, as shown in Fig.5.5(a) for
$T=0.8$. Again, this is related to the local moment like spectra at small
$S\lesssim 0.2$ [Figs.A.1(a), A.2(a)]. On the other hand, $v_{B}$ increases
exponentially with $S$ in the SG phase as shown in Fig.5.5(b), in contrast to
power-law decay of $\lambda_{\mathrm{L}}(S)$ [ig.4.4(b)]. We also plot the
proxy for the diffusion constant $\tilde{D}\sim
v_{B}^{2}/\lambda_{\mathrm{L}}$ as a function of $S$ in Fig.5.6. $\tilde{D}$
increases monotonically with $S$ from the PM to SG phase across the phase
transition.

Figure 5.8: Temperature dependence of $v_{B}$ in the SG and PM phases of SY
model. (a) $v_{B}(S,T)$ (circles) in the SG phase is fitted with
$v_{B}(S,T)\simeq v_{B}(S,0)+aT^{\beta_{v}}$ (line) for $S=1.0$ where
$v_{B}(S,0)\simeq 0.46$ and $\beta_{v}\simeq 1.6$. (b) Logarithmic fit (line)
to $v_{B}(T)$ in the PM phase, $v_{B}\sim a\ln{T}+c$ is shown on a semi-log
plot ($v_{B}$ vs. $\ln{T}$) for $S=1.0$. Figure 5.9: $\tilde{D}\sim
v_{B}^{2}/\lambda_{\mathrm{L}}$ as a function of temperature for $S=0.5$ and
$S=1.0$ in SY model. The SG transition temperatures are shown by vertical
dashed lines.

Temperature dependence of butterfly velocity.— $v_{B}$ is shown as a function
of $T$ for $S=1.0$ across the PM-SG transition. We show $v_{B}(T)$ for $S=0.5$
only in the PM phase, as we did compute $v_{B}$ for the narrow temperature
range $T<T_{SG}$ in the SG phase for $S=0.5$. The butterfly velocity
monotonically increases with temperature through $T_{SG}$, saturating at high
temperature, similar to $\lambda_{\mathrm{L}}(T)$ [Fig.4.5]. As seen in
Fig.5.7, $v_{B}$ approaches a non-zero value for $T\to 0$ in the SG phase.
This implies a finite speed $v_{B}(S,0)$ of information scrambling at $T=0$ in
the SG phase of the SY chain. As shown in Fig.5.8(a), $v_{B}(T)-v_{B}(S,0)\sim
T^{1.6}$, i.e., $v_{B}$ increases with temperature as a power law from the
$T=0$ value. Like the $\lambda_{\mathrm{L}}(T)$ [Fig.5.8(b)], the $T$
dependence of $v_{B}$ is logarithmic, $\sim a\ln{T}+c$, for $T\gtrsim T_{SG}$
in the PM phase [Fig.5.8(b)]. The proxy diffusion constant $\tilde{D}\sim
v_{B}^{2}/\lambda_{\mathrm{L}}$ is plotted in Fig.5.9. For $S=1.0$,
$\tilde{D}$ approaches a large value at low temperature $T\ll T_{SG}$ in the
SG phase, since $v(S,T)\to v_{B}(S,0)\neq 0$ and $\lambda_{\mathrm{L}}\sim T$
for $T\to 0$. In the SG phase one naively expects the transport to slow down
and $\tilde{D}$ to decrease with temperature. Thus, it appears that a system
made out of spin glass quantum dots, like in the SY chain [Eq.(2.1)], can lead
to unusual enhanced diffusion at low temperature. However, this needs to be
verified by calculating actual diffusion coefficient $D$, instead of
$\tilde{D}$. The latter also has dip at the transition and saturates at high
temperature in the PM phase.

Figure 5.10:  Diagrammatic representation of the Kernel equation in Eq:5.13
for $\mathcal{O}(1/N)$ term, $\mathcal{F}_{x}^{a}$ in the OTOC
$F_{x}^{a}(t_{1},t_{2})$. The solid horizontal lines denote dressed retarded
and advanced propagators $G_{R}$, $G_{A}$ and the vertical rung denotes the
Wightmann correlation $G_{lr}$ (Appendix B). The dashed and dotted lines
denote disorder averaging over $J$ and $J^{+}$ couplings in Eq.(2.2). Here
$x^{\prime}=x\pm 1$.

###  5.2 Butterfly velocity across spin glass transition in the $p$-spin glass
chain

To compute the $v_{B}(\Gamma,T)$ in the $p$-spin glass chain [Eq.(2.2)], we
use methods similar to that for the SY chain discussed above. We calculate the
following regularized OTOC defined as,

$\displaystyle F_{x}(t_{1},t_{2})=$
$\displaystyle\frac{1}{N^{2}}\sum_{i,j}\overline{\mathrm{Tr}\Big{[}ys_{i,x}(t_{1})ys_{j,0}(0)ys_{i,x}(t_{2})ys_{j,0}(0)\Big{]}}$
( 5.12)

Following the calculations discussed earlier for the SY chain, the diagrams
that contribute to the OTOC can be grouped in powers of $1/N$. The leading
order $\mathcal{O}(1)$ diagrams are disconnected and do not contribute to the
growth. At the next leading order $\mathcal{O}(1/N)$ we have “ladder” diagrams
which contribute to the growth of OTOC in the scrambling regime,
$\lambda_{\mathrm{L}}^{-1}\lesssim t\lesssim\lambda_{\mathrm{L}}^{-1}\ln{N}$.
As earlier, these ladder diagrams can be written in the form of a Bethe-
Salpeter equation shown in Fig.5.10, i.e.

$\displaystyle\mathcal{F}^{a}_{x}(t_{1},t_{2})=$ $\displaystyle\int
dx^{\prime}\int
dt_{3}dt_{4}K^{a}_{xx^{\prime}}(t_{1},t_{2},t_{3},t_{4})\mathcal{F}^{a}_{x^{\prime}}(t_{3},t_{4}),$
( 5.13)

where the kernel is given by

$\displaystyle K^{a}_{xx^{\prime}}(t_{1},t_{2},t_{3},t_{4})=$
$\displaystyle\left[(3J^{2}+{J^{\prime}}^{2})\delta_{x^{\prime}x}+2{J^{\prime}}^{2}\delta_{x^{\prime},x\pm
1}\right]$ $\displaystyle G_{R}(t_{13})G_{R}(t_{24})G_{W}(t_{34})$ ( 5.14)

After Fourier transforming $\mathcal{F}^{a}_{x}(t_{1},t_{2})$ to momentum
space, we obtain $\mathcal{F}^{a}_{p}(t_{1},t_{2})$ and the Kernel

$\displaystyle
K^{a}_{p}(t_{1},t_{2},t_{3},t_{4})=\tilde{J}^{2}[3+2\alpha(2\cos{p}-1)]$ (
5.15)

where we have used, $J^{2}+{J^{\prime}}^{2}=\tilde{J}^{2}$ and
$\alpha={J^{\prime}}^{2}/\tilde{J}^{2}$.

To calculate the butterfly velocity we again use the real and imaginary $p$
methods discussed for the SY chain above [Sec.5.1]. In the imaginary $p$
method, we employ the single mode ansatz

$\displaystyle
F_{p}^{a}(t_{1},t_{2})\approx\frac{e^{\lambda_{\mathrm{L}}(p)(t_{1}+t_{2})/2}}{C(p)}\Upsilon_{p}^{R}(t_{12})\Upsilon_{p}^{A}(0)$
( 5.16)

We skip the details of the implementations of the two methods here, since they
are similar to those for the SY chain. $v_{B}$ is calculated across the spin
glass phase transition line as a function of $T$ and $\Gamma$.

Figure 5.11: Variation of butterfly velocity with quantum parameter $\Gamma$
in the $p$-spin glass model. $v_{B}$ across the $p$-spin glass phase
transition as a function of $\Gamma$ at $T=0.4$ , $\alpha=0.25$. In the real
$p$ calculation $v_{B}$ is calculated from the slope of light cone for
$\mathcal{F}_{th}=0.5$, for a lattice of $L=50$ sites (see Sec.5.1). The error
bars in the real $p$ results indicate variation of $v_{B}$ over a range of
threshold values $\mathcal{F}_{th}=0.1-1.0$.

$\Gamma$ dependence of butterfly velocity.— Fig.5.11 shows the butterfly
velocity as a function of the quantum fluctuation parameter $\Gamma$ for
$T=0.4$ (in units of $\tilde{J}=1$) ($y$ cut in Fig.3.2). Unlike the SY model
[Fig.5.4], $v_{B}$ does not have a dip at the transition as function of the
quantum parameter. In the PM phase, $v_{B}$ exhibits non-monotonic dependence
on $\Gamma$ with the broad maximum above the glass transition at
$\Gamma=\Gamma_{SG}$, as in the case of $\lambda_{\mathrm{L}}(\Gamma)$
[Fig.4.8].

Figure 5.12: $\Gamma$ dependence of butterfly velocity in the SG and PM phases
of $p$-spin glass model. $v_{B}(\Gamma,T)$ (circles) is fitted (a) in the SG
phase with $v_{B}\simeq v_{B}(0,T)+a\Gamma^{\alpha_{v}}$ (line), where
$\alpha_{v}\simeq 0.5$ and $v_{B}(0,T)\simeq 0.005$ for $T=0.4$, and (b) in
the PM phase with $v_{B}\sim\Gamma^{0.9}$ (line), shown on a log-log plot for
$T=0.4$. Figure 5.13:  $\tilde{D}\sim v_{B}^{2}/\lambda_{\mathrm{L}}$ across
the SG phase transition as a function of $\Gamma$ calculated for $T=0.4$

The butterfly velocity follows a power law,
$v_{B}(\Gamma,T)-v_{B}(0,T)\sim\sqrt{\Gamma}$, with the quantum fluctuation
parameter $\Gamma$ for small $\Gamma<\Gamma_{SG}$ with extrapolated
$v_{B}(0,T)\approx 0$ in the $\Gamma\to 0$ limit, deep in the SG phase
[Fig.5.12(a)]. This is unlike the exponential dependence of $v_{B}$ on
quantum/classical parameter $S$ for the SY chain [Fig.4.4(b)]. $v_{B}$ varies
almost linearly with $\Gamma$ [Fig.5.12(b)] in the PM phase close to the
transition ($\Gamma\gtrsim\Gamma_{SG}$) before reaching the maximum in
Fig.5.11. We also plot the proxy for diffusion constant $\tilde{D}\sim
v_{B}^{2}/\lambda_{\mathrm{L}}$ in Fig.5.13 for $T=0.4$. $\tilde{D}$ is very
small for small $\Gamma$ and then rises to a peak around the transition
$\Gamma_{SG}$, with further non-monotonic dependence on $\Gamma$ in the PM
phase for large $\Gamma>\Gamma_{SG}$.

Figure 5.14: Temperature dependence of butterfly velocity in $p$-spin model.
$v_{B}$ across the phase transition as a function of temperature at
$\Gamma=2.0$, $\alpha=0.25$. The imaginary momentum calculation uses the
single mode ansatz. In the real momentum calculation $v_{B}$ is calculated
from the edge of the light cone for values of $\mathcal{F}_{th}=0.5$, in a
lattice of $L=50$ sites. Variation in $v_{B}$ from nearby values of
$\mathcal{F}_{th}=0.1$ and $\mathcal{F}_{th}=1.0$ is shown as error bar.

Temperature dependence of butterfly velocity.— We demonstrate the overall
temperature dependence of the butterfly velocity going from the SG to the PM
phase through the transition at $T_{SG}$ for $\Gamma=2.0$ in Fig.5.14.
$v_{B}(T)$ exhibits a broad peak in the PM phase, just like
$\lambda_{\mathrm{L}}(T)$ in Fig.4.9. Figs.5.15(a),(b) show that the butterfly
velocity follows $v_{B}(\Gamma,T)-v_{B}(\Gamma,0)\sim T^{1.5}$, with
substantial extrapolated zero-temperature value $v_{B}(\Gamma,0)$ in the SG
phase, and logarithmic $T$ dependence in the PM phase for $T\gtrsim T_{SG}$,
as in the SY model [Fig.5.8]. The variation of resultant $\tilde{D}\sim
v_{B}^{2}/\lambda_{\mathrm{L}}$ with temperature [Fig.5.16] is also very
similar to that for the SY model, except, for the former there is no dip of
$\tilde{D}$ at the transition, unlike the dip in $\tilde{D}$ for SY model,
shown in Fig.5.9.

Figure 5.15: Temperature dependence of $v_{B}$ in PM and SG phases of $p$-spin
glass model.. (a) $v_{B}$ vs. $T$ (circles) in SG phase is fitted as
$v_{B}\simeq v_{B}(\Gamma,0)+aT^{\beta_{v}}$ (line) for $\Gamma=2.0$, where
$\beta_{v}\simeq 1.45$ and $v_{B}(\Gamma,0)\simeq 0.43$. (b) Logarithmic
dependence of $v_{B}$ on $T$ (circles) in PM phase, shown as a straight line
fit on a semi-log plot for $\Gamma=2.0$. Figure 5.16: $\tilde{D}\sim
v_{B}^{2}/\lambda_{\mathrm{L}}$ across the phase transition as a function of
$T$ calculated for $\Gamma=2.0$. in $p$-spin glass model

##  6 Discussion

In this work, we have studied how quantum and thermal fluctuations influence
propagation of quantum information, as diagnosed by butterfly velocity, along
with chaotic growth rate characterized by a Lyapunov exponent, through one-
dimensional generalizations of two paradigmatic zero-dimensional quantum spin
glass models, Sachdev-Ye and $p$-spin glass models. These lattice models allow
us to map out variations of $\lambda_{\mathrm{L}}$ and $v_{B}$ with
temperature and quantum parameter deep inside the phases and near the non-
trivial replica symmetry breaking spin glass transitions, as well as through
crossovers between competing paramagnetic phases, e.g., spin liquid and local
moment behaviours, and across the SG transition. We find distinct evolution of
glassy complexity, indicated by maximum, or lack thereof, in the Lyapunov
exponent and butterfly velocity in the PM phase, as a function quantum and
thermal fluctuations in the SY model. The maximum is seen only as a function
of quantum parameter $S$, around the onset of two-step glassy relaxation, away
from the glass transition from the PM side. The maximum is absent when the
glass transition is approached by decreasing $T$. In contrast, the maximum
occurs as function of both $T$ and the quantum parameter $\Gamma$ in the PM
phase of $p$-spin glass. The maximum corresponds to intensified chaos due to
maximal complexity arising from the sampling of exponentially large number of
saddle points of the underlying glassy landscape. In the semiclassical limit,
the maximal number of saddle points leads to extreme sensitivity to the
initial condition Correale  _et al._ (2023) or chaos.

Our works needs to be extended in future in several directions. To better
understand connection between glassy complexity Correale  _et al._ (2023);
Anous and Haehl (2021) and quantum many-body chaos, it will be worthwhile to
compute the complexity for the SY model using both semiclassical methods
Correale  _et al._ (2023) and quantum Thouless-Anderson-Palmer (TAP) equations
Biroli and Cugliandolo (2001), as have been done for quantum $p$-spin glass.
The asymptotic temperature dependence of $\lambda_{\mathrm{L}}$ and $v_{B}$ we
obtain deep in the SG phase of SY and $p$-spin glass implies a finite $v_{B}$
for the $T=0$ SG ground state, and, consequently, lead to rapid increase of
the quantity $\tilde{D}\sim v_{B}^{2}/\lambda_{\mathrm{L}}$ for $T\to 0$. As
$\tilde{D}$ often acts as a proxy to diffusion constant, it will be
interesting to probe the possibility of such anomalously large diffusion
through actual transport calculations in a chain of spin glass quantum dots.

Our real-frequency saddle-point calculations of $\lambda_{\mathrm{L}}$ and
$v_{B}$, and corresponding bosonic spectral functions, in the SY model
indicate subtle competition between spin liquid and local moment solutions, as
well as strong $S$ and $T$ dependent corrections to the conformal solution for
the SYK-like spin liquid maximal chaotic behaviour, even at small values of
$S\lesssim 0.05$ and low temperature $T\sim 0.1$. This implies intriguing
crossover in the small $S$ region towards the expected maximally chaotic SYK
spin liquid for $S\to 0$. We have not studied this crossover and the relative
stability and intricate competition between spin liquid and local moment
solutions in the small $S$ region in much detail in this work. It will be
important to study the chaos and spectral properties in the small-$S$ regime
in future work. Finally, it will be an exciting research direction to compute
entanglement properties in various crossover regions and across the replica
symmetry breaking spin glass transition in these large-$N$ models following
methods Haldar  _et al._ (2020); Bera  _et al._ (2023) similar to those
applied for calculating entanglement entropy in lattice of SYK dots and in
Hubbard model within dynamical mean field theory.

## Acknowldegements

VL acknowledges support from CSIR, India and is grateful to the workstation
provided by CHEP on which most of calculations in this work are done. SB
acknowledges support from SERB (Grant No. CRG/2022/001062), DST, India, and
QuST project of the DST, Govt. of India.

## References

  * Sekino and Susskind (2008) Yasuhiro Sekino and L. Susskind, “Fast scramblers,” Journal of High Energy Physics 2008, 065 (2008), arXiv:0808.2096 [hep-th] .
  * Maldacena  _et al._ (2016) Juan Maldacena, Stephen H. Shenker, and Douglas Stanford, “A bound on chaos,” Journal of High Energy Physics 2016, 1–17 (2016).
  * (3) A. Kitaev, “A simple model of quantum holography,” Talks at KITP, April 7, 2015 and May 27, 2015 .
  * Sahu  _et al._ (2019) Subhayan Sahu, Shenglong Xu, and Brian Swingle, “Scrambling dynamics across a thermalization-localization quantum phase transition,” Phys. Rev. Lett. 123, 165902 (2019).
  * Jensen (2016) Kristan Jensen, “Chaos in ${\mathrm{ads}}_{2}$ holography,” Phys. Rev. Lett. 117, 111601 (2016).
  * Maldacena and Stanford (2016) Juan Maldacena and Douglas Stanford, “Remarks on the sachdev-ye-kitaev model,” Phys. Rev. D 94, 106002 (2016).
  * Patel and Sachdev (2017) Aavishkar A. Patel and Subir Sachdev, “Quantum chaos on a critical fermi surface,” Proceedings of the National Academy of Sciences 114, 1844–1849 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1618185114 .
  * Banerjee and Altman (2017) Sumilan Banerjee and Ehud Altman, “Solvable model for a dynamical quantum phase transition from fast to slow scrambling,” Phys. Rev. B 95, 134302 (2017).
  * Chowdhury  _et al._ (2022) Debanjan Chowdhury, Antoine Georges, Olivier Parcollet, and Subir Sachdev, “Sachdev-ye-kitaev models and beyond: Window into non-fermi liquids,” Rev. Mod. Phys. 94, 035004 (2022).
  * Srednicki (1994) Mark Srednicki, “Chaos and quantum thermalization,” Phys. Rev. E 50, 888–901 (1994).
  * Murthy and Srednicki (2019) Chaitanya Murthy and Mark Srednicki, “Bounds on chaos from the eigenstate thermalization hypothesis,” Phys. Rev. Lett. 123, 230606 (2019).
  * Aleiner  _et al._ (2016) Igor L. Aleiner, Lara Faoro, and Lev B. Ioffe, “Microscopic model of quantum butterfly effect: Out-of-time-order correlators and traveling combustion waves,” Annals of Physics 375, 378–406 (2016).
  * Stanford (2016) Douglas Stanford, “Many-body chaos at weak coupling,” Journal of High Energy Physics 2016, 9 (2016).
  * Rozenbaum  _et al._ (2017) Efim B. Rozenbaum, Sriram Ganeshan, and Victor Galitski, “Lyapunov exponent and out-of-time-ordered correlator’s growth rate in a chaotic system,” Phys. Rev. Lett. 118, 086801 (2017).
  * Scaffidi and Altman (2019) Thomas Scaffidi and Ehud Altman, “Chaos in a classical limit of the Sachdev-Ye-Kitaev model,” Phys. Rev. B 100, 155128 (2019).
  * Xu  _et al._ (2020) Tianrui Xu, Thomas Scaffidi, and Xiangyu Cao, “Does scrambling equal chaos?” Phys. Rev. Lett. 124, 140602 (2020).
  * Kob (1997) W. Kob, “The Mode-Coupling Theory of the Glass Transition,” arXiv e-prints , cond-mat/9702073 (1997), arXiv:cond-mat/9702073 [cond-mat.stat-mech] .
  * Debenedetti and Stillinger (2001) Pablo G. Debenedetti and Frank H. Stillinger, “Supercooled liquids and the glass transition,” Nature 410, 259–267 (2001).
  * Berthier and Biroli (2011) Ludovic Berthier and Giulio Biroli, “Theoretical perspective on the glass transition and amorphous materials,” Rev. Mod. Phys. 83, 587–645 (2011).
  * Büchner and Heuer (1999) Stephan Büchner and Andreas Heuer, “Potential energy landscape of a model glass former: Thermodynamics, anharmonicities, and finite size effects,” Phys. Rev. E 60, 6507–6518 (1999).
  * Sastry (2001) Srikanth Sastry, “The relationship between fragility, configurational entropy and the potential energy landscape of glass-forming liquids,” Nature (London) 409, 164–167 (2001), arXiv:cond-mat/0011317 [cond-mat.stat-mech] .
  * Sciortino (2005) Francesco Sciortino, “Potential energy landscape description of supercooled liquids and glasses,” Journal of Statistical Mechanics: Theory and Experiment 2005, P05015 (2005).
  * Binder and Young (1986) K. Binder and A. P. Young, “Spin glasses: Experimental facts, theoretical concepts, and open questions,” Rev. Mod. Phys. 58, 801–976 (1986).
  * Vincent (2007) E. Vincent, “Ageing, rejuvenation and memory: The example of spin-glasses,” in  _Ageing and the Glass Transition_, edited by Malte Henkel, Michel Pleimling, and Roland Sanctuary (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 7–60.
  * Hodgdon and Stillinger (1993) Jennifer A. Hodgdon and Frank H. Stillinger, “Stokes-einstein violation in glass-forming liquids,” Phys. Rev. E 48, 207–213 (1993).
  * Karmakar  _et al._ (2010) Smarajit Karmakar, Chandan Dasgupta, and Srikanth Sastry, “Analysis of dynamic heterogeneity in a glass former from the spatial correlations of mobility,” Phys. Rev. Lett. 105, 015701 (2010).
  * Edwards and Anderson (1975) S. F. Edwards and P. W. Anderson, “Theory of spin glasses,” Journal of Physics F Metal Physics 5, 965–974 (1975).
  * Mezard  _et al._ (1987) Marc Mezard, Giorgio Parisi, and Miguel Angel Virasoro,  _Spin Glass Theory and Beyond: AN Introduction to the Replica Method and its Applications_ (1987).
  * Parisi (1979) G. Parisi, “Infinite number of order parameters for spin-glasses,” Phys. Rev. Lett. 43, 1754–1756 (1979).
  * Castellani and Cavagna (2005) Tommaso Castellani and Andrea Cavagna, “Spin-glass theory for pedestrians,” Journal of Statistical Mechanics: Theory and Experiment 2005, P05012 (2005).
  * Bera  _et al._ (2022) Surajit Bera, K. Y. Venkata Lokesh, and Sumilan Banerjee, “Quantum-to-classical crossover in many-body chaos and scrambling from relaxation in a glass,” Phys. Rev. Lett. 128, 115302 (2022).
  * Nieuwenhuizen (1995) Th. M. Nieuwenhuizen, “Exactly solvable model of a quantum spin glass,” Phys. Rev. Lett. 74, 4289–4292 (1995).
  * Nieuwenhuizen and Ritort (1998) Theo M. Nieuwenhuizen and Felix Ritort, “Quantum phase transition in spin glasses with multi-spin interactions,” Physica A Statistical Mechanics and its Applications 250, 8–45 (1998), arXiv:cond-mat/9706244 [cond-mat.dis-nn] .
  * Cugliandolo  _et al._ (2000) Leticia F. Cugliandolo, Daniel R. Grempel, and Constantino A. da Silva Santos, “From second to first order transitions in a disordered quantum magnet,” Phys. Rev. Lett. 85, 2589–2592 (2000).
  * Correale  _et al._ (2023) Lorenzo Correale, Anatoli Polkovnikov, Marco Schirò, and Alessandro Silva, “Probing chaos in the spherical p-spin glass model,” arXiv e-prints , arXiv:2303.15393 (2023), arXiv:2303.15393 [cond-mat.dis-nn] .
  * Auffinger and Ben Arous (2011) Antonio Auffinger and Gerard Ben Arous, “Complexity of random smooth functions on the high-dimensional sphere,” arXiv e-prints , arXiv:1110.5872 (2011), arXiv:1110.5872 [math.PR] .
  * Anous and Haehl (2021) Tarek Anous and Felix M Haehl, “The quantum p-spin glass model: a user manual for holographers,” Journal of Statistical Mechanics: Theory and Experiment 2021, 113101 (2021).
  * Sachdev and Ye (1993) Subir Sachdev and Jinwu Ye, “Gapless spin-fluid ground state in a random quantum heisenberg magnet,” Phys. Rev. Lett. 70, 3339–3342 (1993).
  * Georges  _et al._ (2000) Antoine Georges, Olivier Parcollet, and Subir Sachdev, “Mean field theory of a quantum heisenberg spin glass,” Phys. Rev. Lett. 85, 840–843 (2000).
  * Georges  _et al._ (2001) A. Georges, O. Parcollet, and S. Sachdev, “Quantum fluctuations of a nearly critical heisenberg spin glass,” Phys. Rev. B 63, 134406 (2001).
  * Gu  _et al._ (2017) Yingfei Gu, Xiao-Liang Qi, and Douglas Stanford, “Local criticality, diffusion and chaos in generalized sachdev-ye-kitaev models,” Journal of High Energy Physics 2017, 125 (2017).
  * Sachdev (2015) Subir Sachdev, “Bekenstein-hawking entropy and strange metals,” Phys. Rev. X 5, 041025 (2015).
  * Song  _et al._ (2017a) Xue-Yang Song, Chao-Ming Jian, and Leon Balents, “Strongly correlated metal built from sachdev-ye-kitaev models,” Phys. Rev. Lett. 119, 216601 (2017a).
  * Song  _et al._ (2017b) Xue-Yang Song, Chao-Ming Jian, and Leon Balents, “Strongly correlated metal built from sachdev-ye-kitaev models,” Phys. Rev. Lett. 119, 216601 (2017b).
  * Davison  _et al._ (2017) Richard A. Davison, Wenbo Fu, Antoine Georges, Yingfei Gu, Kristan Jensen, and Subir Sachdev, “Thermoelectric transport in disordered metals without quasiparticles: The sachdev-ye-kitaev models and holography,” Phys. Rev. B 95, 155131 (2017).
  * Jian  _et al._ (2017) Chao-Ming Jian, Zhen Bi, and Cenke Xu, “Model for continuous thermal metal to insulator transition,” Phys. Rev. B 96, 115122 (2017).
  * Jian and Yao (2017) Shao-Kai Jian and Hong Yao, “Solvable sachdev-ye-kitaev models in higher dimensions: From diffusion to many-body localization,” Phys. Rev. Lett. 119, 206602 (2017).
  * Berkooz  _et al._ (2017) Micha Berkooz, Prithvi Narayan, Moshe Rozali, and Joan Simón, “Higher dimensional generalizations of the SYK model,” Journal of High Energy Physics 2017, 138 (2017), arXiv:1610.02422 [hep-th] .
  * Zhang (2017) Pengfei Zhang, “Dispersive sachdev-ye-kitaev model: Band structure and quantum chaos,” Phys. Rev. B 96, 205138 (2017).
  * Chowdhury  _et al._ (2018) Debanjan Chowdhury, Yochai Werman, Erez Berg, and T. Senthil, “Translationally invariant non-fermi-liquid metals with critical fermi surfaces: Solvable models,” Phys. Rev. X 8, 031024 (2018).
  * Haldar  _et al._ (2018) Arijit Haldar, Sumilan Banerjee, and Vijay B. Shenoy, “Higher-dimensional sachdev-ye-kitaev non-fermi liquids at lifshitz transitions,” Phys. Rev. B 97, 241106 (2018).
  * Guo  _et al._ (2019) Haoyu Guo, Yingfei Gu, and Subir Sachdev, “Transport and chaos in lattice sachdev-ye-kitaev models,” Phys. Rev. B 100, 045140 (2019).
  * Chowdhury and Swingle (2017) Debanjan Chowdhury and Brian Swingle, “Onset of many-body chaos in the $o(n)$ model,” Phys. Rev. D 96, 065005 (2017).
  * Sahu and Swingle (2020) Subhayan Sahu and Brian Swingle, “Information scrambling at finite temperature in local quantum systems,” Phys. Rev. B 102, 184303 (2020).
  * Patel  _et al._ (2017) Aavishkar A. Patel, Debanjan Chowdhury, Subir Sachdev, and Brian Swingle, “Quantum butterfly effect in weakly interacting diffusive metals,” Phys. Rev. X 7, 031047 (2017).
  * Bilitewski  _et al._ (2018) Thomas Bilitewski, Subhro Bhattacharjee, and Roderich Moessner, “Temperature dependence of the butterfly effect in a classical many-body system,” Phys. Rev. Lett. 121, 250602 (2018).
  * Ruidas and Banerjee (2021) Sibaram Ruidas and Sumilan Banerjee, “Many-body chaos and anomalous diffusion across thermal phase transitions in two dimensions,” SciPost Physics 11, 087 (2021), arXiv:2007.12708 [cond-mat.stat-mech] .
  * Bilitewski  _et al._ (2021) Thomas Bilitewski, Subhro Bhattacharjee, and Roderich Moessner, “Classical many-body chaos with and without quasiparticles,” Phys. Rev. B 103, 174302 (2021).
  * Blake (2016) Mike Blake, “Universal diffusion in incoherent black holes,” Phys. Rev. D 94, 086014 (2016).
  * Swingle and Chowdhury (2017) Brian Swingle and Debanjan Chowdhury, “Slow scrambling in disordered quantum systems,” Phys. Rev. B 95, 060201 (2017).
  * Wei  _et al._ (2019) Bo-Bo Wei, Gaoyong Sun, and Myung-Joong Hwang, “Dynamical scaling laws of out-of-time-ordered correlators,” Phys. Rev. B 100, 195107 (2019).
  * Gu and Kitaev (2019) Yingfei Gu and Alexei Kitaev, “On the relation between the magnitude and exponent of OTOCs,” JHEP 02, 075 (2019), arXiv:1812.00120 [hep-th] .
  * Biroli and Parcollet (2002) Giulio Biroli and Olivier Parcollet, “Out-of-equilibrium dynamics of a quantum heisenberg spin glass,” Phys. Rev. B 65, 094414 (2002).
  * Grempel and Rozenberg (1998) D. R. Grempel and M. J. Rozenberg, “Fluctuations in a quantum random heisenberg paramagnet,” Phys. Rev. Lett. 80, 389–392 (1998).
  * Cugliandolo  _et al._ (2001) Leticia F. Cugliandolo, D. R. Grempel, and Constantino A. da Silva Santos, “Imaginary-time replica formalism study of a quantum spherical p-spin-glass model,” Phys. Rev. B 64, 014403 (2001).
  * Götze (2008) W. Götze,  _Complex Dynamics of Glass-Forming Liquids_ (Oxford University Press, Oxford, 2008).
  * Reichman and Charbonneau (2005) David R. Reichman and Patrick Charbonneau, “Mode-coupling theory,” Journal of Statistical Mechanics: Theory and Experiment 2005, 05013 (2005), arXiv:cond-mat/0511407 [cond-mat.soft] .
  * Camjayi and Rozenberg (2003) Alberto Camjayi and Marcelo J. Rozenberg, “Quantum and thermal fluctuations in the su($n$) heisenberg spin-glass model near the quantum critical point,” Phys. Rev. Lett. 90, 217202 (2003).
  * Christos  _et al._ (2022) Maine Christos, Felix M. Haehl, and Subir Sachdev, “Spin liquid to spin glass crossover in the random quantum heisenberg magnet,” Phys. Rev. B 105, 085120 (2022).
  * Arrachea and Rozenberg (2002) Liliana Arrachea and Marcelo J. Rozenberg, “Infinite-range quantum random heisenberg magnet,” Phys. Rev. B 65, 224430 (2002).
  * Cugliandolo and Lozano (1999) Leticia F. Cugliandolo and Gustavo Lozano, “Real-time nonequilibrium dynamics of quantum glassy systems,” Phys. Rev. B 59, 915–942 (1999).
  * Gu  _et al._ (2020) Yingfei Gu, Alexei Kitaev, Subir Sachdev, and Grigory Tarnopolsky, “Notes on the complex Sachdev-Ye-Kitaev model,” JHEP 02, 157 (2020), arXiv:1910.14099 [hep-th] .
  * Kurchan (2016) Jorge Kurchan, “Quantum bound to chaos and the semiclassical limit,” arXiv e-prints , arXiv:1612.01278 (2016), arXiv:1612.01278 [cond-mat.stat-mech] .
  * Khemani  _et al._ (2018) Vedika Khemani, David A. Huse, and Adam Nahum, “Velocity-dependent lyapunov exponents in many-body quantum, semiclassical, and classical chaos,” Phys. Rev. B 98, 144304 (2018).
  * Roberts and Swingle (2016) Daniel A. Roberts and Brian Swingle, “Lieb-robinson bound and the butterfly effect in quantum field theories,” Phys. Rev. Lett. 117, 091602 (2016).
  * Lieb and Robinson (1972) Elliott H. Lieb and Derek W. Robinson, “The finite group velocity of quantum spin systems,” Communications in Mathematical Physics 28, 251–257 (1972).
  * Kitaev and Suh (2018) Alexei Kitaev and S. Josephine Suh, “The soft mode in the Sachdev-Ye-Kitaev model and its gravity dual,” JHEP 05, 183 (2018), arXiv:1711.08467 [hep-th] .
  * Biroli and Cugliandolo (2001) Giulio Biroli and Leticia F. Cugliandolo, “Quantum thouless-anderson-palmer equations for glassy systems,” Phys. Rev. B 64, 014206 (2001).
  * Haldar  _et al._ (2020) Arijit Haldar, Surajit Bera, and Sumilan Banerjee, “Rényi entanglement entropy of fermi and non-fermi liquids: Sachdev-ye-kitaev model and dynamical mean field theories,” Phys. Rev. Res. 2, 033505 (2020).
  * Bera  _et al._ (2023) Surajit Bera, Arijit Haldar, and Sumilan Banerjee, “Dynamical mean-field theory for Rényi entanglement entropy and mutual Information in Hubbard Model,” arXiv e-prints , arXiv:2302.10940 (2023), arXiv:2302.10940 [cond-mat.str-el] .
  * Parcollet  _et al._ (1998) Olivier Parcollet, Antoine Georges, Gabriel Kotliar, and Anirvan Sengupta, “Overscreened multichannel $\mathrm{SU}(n)$ kondo model: Large-$n$ solution and conformal field theory,” Phys. Rev. B 58, 3794–3813 (1998).

##  Appendix A Imaginary-time path integral and saddle-point equations for SY
chain

We write the partition function
$\mathcal{Z}(\\{J_{ij,x},J^{\prime}_{ij,x}\\})=\mathrm{Tr}[\exp(-\beta\mathcal{H})]$
($\beta=1/T$) of the disordered model of Eq.(2.1) as an imaginary-time path
integral using the bosonic representation $S_{i\alpha\beta x}=b_{i\alpha
x}^{\dagger}b_{i\beta x}-S\delta_{\alpha\beta}$ with the constraint
$\sum_{\alpha}b_{i\alpha x}^{\dagger}b_{i\alpha x}=SM$ ($S\geq 0$) on the
number of bosons at each $i,x$, to fix the spin of the bosonic representation
to $S$. Here the (real) random couplings $\\{J_{ij,x},J^{\prime}_{ij,x}\\}$
are drawn from Gaussian distributions with zero mean and variances $J^{2}$ and
$J^{\prime 2}$, respectively. We calculate disorder-averaged quantities by
introducing replicas $a=1,\dots,n$ to obtain
$\overline{\mathcal{Z}^{n}}=\int\mathcal{D}(\bar{b},b)\mathcal{D}\lambda\mathrm{exp}[-S_{eff}]$,
where

$\displaystyle S_{eff}=$
$\displaystyle\int_{0}^{\beta}d\tau\sum_{ix}\left[\bar{b}_{i\alpha
x}^{a}\partial_{\tau}b_{i\alpha x}^{a}-\lambda^{a}_{ix}(\tau)(\bar{b}_{i\alpha
x}^{a}b_{i\alpha x}^{a}-SM)\right]$ $\displaystyle-\frac{1}{2MN}$
$\displaystyle\int_{0}^{\beta}d\tau