# The Lyapunov exponent as a signature of dissipative
many-body quantum chaos

Antonio M. Garc\́mathrm{i}a-Garc\́mathrm{i}a<EMAIL_ADDRESS>Shanghai Center
for Complex Physics, School of Physics and Astronomy, Shanghai Jiao Tong
University, Shanghai 200240, China  Jacobus J. M. Verbaarschot
<EMAIL_ADDRESS>Center for Nuclear Theory and Department
of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11794,
USA  Jie-ping Zheng (郑杰平<EMAIL_ADDRESS>cn Shanghai Center for Complex
Physics, School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China

###### Abstract

A distinct feature of Hermitian quantum chaotic dynamics is the exponential
increase of certain out-of-time-order-correlation (OTOC) functions around the
Ehrenfest time with a rate given by a Lyapunov exponent. Physically, the OTOCs
describe the growth of quantum uncertainty that crucially depends on the
nature of the quantum motion. Here, we employ the OTOC in order to provide a
precise definition of dissipative quantum chaos. For this purpose, we compute
analytically the Lyapunov exponent for the vectorized formulation of the large
$q$-limit of a $q$-body Sachdev-Ye-Kitaev model coupled to a Markovian bath.
These analytic results are confirmed by an explicit numerical calculation of
the Lyapunov exponent for several values of $q\geq 4$ based on the solutions
of the Schwinger-Dyson and Bethe-Salpeter equations. We show that the Lyapunov
exponent decreases monotonically as the coupling to the bath increases and
eventually becomes negative at a critical value of the coupling signaling a
transition to a dynamics which is no longer quantum chaotic. Therefore, a
positive Lyapunov exponent is a defining feature of dissipative many-body
quantum chaos. The observation of the breaking of the exponential growth for
sufficiently strong coupling suggests that dissipative quantum chaos may
require in certain cases a sufficiently weak coupling to the environment.

###### Contents

  1. I Introduction
  2. II The dissipative SYK model and the calculation of Green’s functions
    1. II.1 SD equations and calculation of Green’s functions
  3. III Calculation of the OTOC
    1. III.1 Numerical OTOC and Lyapunov exponent
    2. III.2 Analytic calculation of $\lambda_{L}$ in the large $q$ limit
  4. IV Discussion and Conclusions
  5. A Conventions
  6. B Solution of the Schwinger-Dyson Equations
  7. C Derivation of the OTOC kernel equation
  8. D Solution of kernel equation for $\beta=0$

##  I Introduction

Research on many-body quantum chaos is receiving a lot of attention because
its strikingly universal features allow a qualitative description of the
quantum time evolution that does not depend on the details of the Hamiltonian.
These universal features leave fingerprints in the quantum dynamics at
different time scales.

One of these time scales is the so called Ehrenfest time, or scrambling [1, 2]
time, $t_{E}$, characterized by the time scale below which the quantum
dynamics largely follows the classical motion. For $t\gtrsim t_{E}$,
deviations from the classical dynamics due to the increase of quantum
uncertainty develop rapidly. For semiclassical quantum chaotic systems, this
is reflected in the exponential growth of certain out-of-time correlation
(OTOC) functions [3, 4] with a rate controlled by the classical Lyapunov
exponent describing the exponential divergence of classical trajectories. By
contrast, for integrable systems, the growth of the OTOC is generically slower
than exponential. This makes the calculation of the OTOC a powerful tool to
characterize the nature of the quantum dynamics around the relatively short
time scale of the Ehrenfest time.

For much longer times, of the order of the Heisenberg time, the inverse of the
mean level spacing, quantum chaotic dynamics is well characterized by level
statistics. According to the celebrated Bohigas-Giannoni-Schmit (BGS)
conjecture [5], level statistics of quantum chaotic systems agrees with the
predictions of random matrix theory [6, 7, 8, 9, 10, 11, 12]. The study of
spectral correlations has been by far the most used tool to study quantum
chaos because it is relatively easy to obtain the spectrum numerically and
compare it with the compact analytic results provided by random matrix theory
[12].

However, in the early days of the development of quantum chaos theory, the
OTOC, and the Lyapunov exponent, were computed analytically, or semi-
analytically, for simple systems such as a kicked rotor [4] or a non-
interacting particle in a random potential [3]. For many-body quantum chaotic
systems, it was considered quite challenging to compute the OTOC either
numerically or analytically.

Recently, the situation changed drastically after two major results: the
proposal that in quantum chaotic systems the Lyapunov exponent is subject to a
universal bound [13] which is conjectured to be saturated in field theories
with a gravity dual, and the observation [14] that, in the low temperature
limit, the bound is saturated in the now called Sachdev-Ye-Kitaev (SYK) model
[15, 16, 17, 18, 19, 20, 21, 14, 22], $N$ zero-dimensional Majorana fermions
with random $q$-body interactions of infinite range in Fock space. Moreover,
recent developments in numerical techniques have made it possible to largely
confirm these theoretical predictions [23, 24].

Therefore, the SYK model is quantum chaotic at the relatively short time scale
of the Ehrenfest time. It is indeed quantum chaotic at all scales as spectral
correlations [25] are well described by random matrix theory.

All the above results apply only to Hermitian systems. However, the dynamics
of a quantum chaotic system coupled to an environment, so that the system
Hamiltonian is not Hermitian, has received a recent boost of attention,
especially for systems that are PT symmetric [26], in a broad variety of
problems from quantum optics [27, 28] and condensed matter [29, 30], to cold
atom [31], and high energy physics [32], see Ref. [33] for a comprehensive
review.

However, the meaning and characterization of quantum chaos in these systems
are much less developed. Although the BGS conjecture has effectively been
extended to non-Hermitian systems [34, 35], and there have been several
successful comparisons between quantum dissipative systems [36, 37, 38, 39,
40, 41, 42] with the corresponding non-Hermitian random matrix ensembles [43,
44, 45], it still remains poorly understood the degree of universality [46] of
the conjecture in non-Hermitian systems and, more generally, the fate of
quantum chaotic effects in the presence of an environment. For instance, one
would expect that a strong coupling to the bath would erase, partially or
completely, quantum features but it is unclear how this would precisely be
reflected in observables like the OTOC.

Here, we address this problem by studying the OTOC of an SYK model coupled to
a Markovian bath [47, 48, 49] in the framework of the so called Lindblad
formalism [50, 51, 52, 53, 54].

We aim to answer the following questions: does the OTOC still grow
exponentially around the Ehrenfest time in the presence of an environment? If
this applies, it could be used as a sharp definition of dissipative quantum
chaos. Does a quantum-classical transition occur for a sufficiently strong
coupling to the environment? If so, this will constrain the definition of
dissipative quantum chaos to the quantum region where the Lyapunov exponent is
positive, no matter whether spectral correlations are described by non-
Hermitian random matrix theory [44].

The role of decoherence in the OTOC [55, 56, 57, 58] and whether chaotic and
integrable motion can be distinguishable in a quantum dissipative system [59]
have already been investigated in the literature but, as far as we know, the
question stated above has not yet been fully settled.

Other studies of OTOCs in non-Hermitian quantum chaotic systems, limited to
single particle kicked rotors [60], two level systems [61], the Ising model
[62], or considering non-Hermitian operators [63] in the OTOC, altogether
avoided addressing the question of the definition of dissipative many-body
quantum chaos that motivates our research. Likewise, different aspects of the
dynamics of non-Hermitian SYK’s [64, 65, 66, 67, 49, 68, 69, 70, 71, 72, 73,
74] that have confirmed agreement with the RMT predictions, have also been the
subject of recent research, but their emphasis was not on the description of
early signatures of dissipative quantum chaos.

The Lyapunov exponent has been computed [75, 76, 77, 78, 79] in several
settings involving two coupled SYK models that, with a certain choice of
operators, could be understood [80] as a single site SYK model coupled to a
bath modeled by the other SYK model. This relation has been explicitly
exploited in Ref. [75] to show that, in agreement with our results, if the SYK
model acting as a bath has a number of Majoranas parametrically larger than
the SYK acting as a system, the Lyapunov exponent vanishes if the coupling to
the bath is sufficiently strong.

The manuscript is organized as follows: we start in section II with the
introduction of the Lindblad SYK model and its main features. All calculations
are performed for the Keldysh formulation of the vectorized form of this
model, which is also discussed in this section. In section III.1, we first
carry out the calculation of the OTOC from the numerical solution of the
Schwinger-Dyson (SD) equations for $q\geq 4$. The resulting Lyapunov exponent
decreases monotonously with the coupling to the bath and eventually vanishes,
indicating the transition to a non quantum-chaotic time evolution. Section
III.2 is devoted to an analytical calculation of the Lyapunov exponent in the
large $q$ limit that is in full agreement with the numerical findings of the
previous section. In view of these results, in section IV we discuss both a
possible definition of quantum chaos in dissipative systems and the
universality of our results. In the appendix, we cover technical aspects of
the calculation of the OTOC such as the derivation of both the Schwinger-Dyson
and the Bethe-Salpeter equations, the latter also termed kernel equation,
necessary to compute the Lyapunov exponent analytically or numerically.

##  II The dissipative SYK model and the calculation of Green’s functions

We study a single-site SYK model coupled to a Markovian environment described
by the Lindblad formalism [48, 47] in which the bath is characterized by jump
operators $L_{i}$ that depend on the single-site SYK Majoranas. The evolution
of the density matrix of this system is dictated by:

$\displaystyle\frac{d\rho}{dt}=\mathcal{L},\qquad\mathcal{L}=-i[H,\rho]+\sum_{i}L_{i}\rho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\rho\\},$ (1)

where $\mathcal{L}$ is termed the Liouvillian. The Hamiltonian $H$ of the SYK
model is given by,

$\displaystyle H=i^{q/2}\sum_{i_{1}<i_{2}<\cdots<i_{q}}J_{i_{1}\cdots
i_{q}}\psi_{i_{1}}\cdots\psi_{i_{q}}$ (2)

with $J_{i_{1}\cdots i_{q}}$ Gaussian random couplings

$\displaystyle\langle J_{i_{1}i_{2}\cdots i_{q}}\rangle=0,\qquad\langle
J_{ijkl}^{2}\rangle=\frac{(q-1)!J^{2}}{N^{q-1}}$ (3)

and $\\{\psi_{i},\psi_{j}\\}=\delta_{ij}$. We consider jump operators,
$L_{i}=\sqrt{\mu}\psi_{i}$, linear in the Majorana fermions. Since
$\psi_{i}^{2}=1/2$, the Liouvillian simplifies to

$\displaystyle\qquad\mathcal{L}=-i[H,\rho]+\mu\sum_{i}\psi_{i}\rho\psi_{i}-\frac{N\mu}{2}\rho.$
(4)

It is clear that an equilibrium solution of the Lindblad equation is given by
the identity matrix

$\displaystyle\rho=\sum_{k}|k\rangle\langle k|,$ (5)

which, for an equilibrium system, is the density matrix at infinite
temperature. Physically, this means that the effect of the coupling to a bath
induces an out of equilibrium dynamics resulting in an increase of temperature
until the final steady state is reached at infinite temperature.

The study of the dynamics of this system requires [81] the vectorization of
the Liouvillian, resulting in a doubling of degrees of freedom, in order to
set up the Keldysh path integral [82]. This doubling is a direct consequence
of the fact that the density matrix in Eq. (1) is an operator that evolves
from the right with the adjoint of the left time evolution, which in the
doubled space can be interpreted as backward time evolution. The Liouvillian
is usually referred to as a super-operator as it generates the dynamics of
another operator, the density matrix. We note that in this formulation, the
infinite temperature thermal density matrix $\rho(\infty)$ is mapped onto the
Thermo-Field-Double (TFD) state at infinite temperature

$\displaystyle|I\rangle\equiv\sum_{k}|k\rangle\otimes|k\rangle.$ (6)

This is an exact eigenstate with a zero eigenvalue describing the steady
($t\to\infty$ limit) state of the vectorized system.

###  II.1 SD equations and calculation of Green’s functions

In order to carry out the mentioned vectorization, following Ref. [49], we
first employ the Choi-Jamiolkowski isomorphism to express the density matrix
operator as a vector in a space where the number of degrees of freedom are
doubled. In a second step, we set up the Keldysh integral by relating
Majoranas in this doubled Fock space to a certain time contour representing
backward and forward time evolution. Vectors in the doubled space,
$\psi\equiv\psi^{+}\otimes\psi^{-}$ obey the following time evolution
equation,

$\displaystyle\frac{d\psi}{dt}=\mathcal{L}\psi,$ (7)

where the vectorized Lindbladian is given by

$\displaystyle\mathcal{L}=-iH^{+}+i(-1)^{\frac{q}{2}}H^{-}-i\mu\sum_{i}\psi^{+}_{i}\psi^{-}_{i}-\frac{1}{2}\mu
N.$ (8)

Here, $H^{+}=H(\psi^{+})\otimes\mathbbm{1}$ and $H^{-}=\mathbbm{1}\otimes
H(\psi^{-})$ are two copies ($\pm$) of the single site SYK Hamiltonian, Eq.
(2), with the same probability distribution but different fermions. For
notational clarity, we are not displaying the tensor products in the remainder
of the paper. All eigenvalues of the Lindblad operator have a non-positive
real part, so that the equilibrium state corresponds to the eigenvalue
$\lambda=0$. In terms of eigenstates of the single SYK Hamiltonian,
$H|k\rangle=E_{k}|k\rangle$, its eigenstate is given by the TFD state Eq. (6).

The path integral for the overlap of the initial state $\psi(t_{i})$ and the
final state $\psi(t_{f})$ is given by

$\displaystyle\langle\psi(t_{f})|\psi(t_{i})\rangle=\int D\psi e^{iS}$ (9)

with action given by

$\displaystyle
iS=\int_{t_{i}}^{t_{f}}dt\left(-\frac{1}{2}\psi_{i}^{-}(t)\partial_{t}\psi_{i}^{-}(t)-\frac{1}{2}\psi_{i}^{+}(t)\partial_{t}\psi_{i}^{+}(t)+\mathcal{L}\right),$
(10)

where we use the notation of summing over repeated indexes. We will be
interested in sufficiently long times where the system is close to a
stationary state. Given the quantum chaotic nature of the dynamics of the SYK
model, we can safely neglect any dependence on initial conditions. Despite the
fact that the system is close to stationary, we still want to evaluate real
time correlation functions for any time difference. As a result of these
considerations, we set $t_{i}=-\infty$, $t_{f}=+\infty$. The Keldysh path
integral formalism involves combining the integrations over $\psi^{+}$ and
$\psi^{-}$ into a single integration variable over a doubled contour in the
complex plane. More specifically, the integrals over $\psi^{+}$ and $\psi^{-}$
are combined into a single contour,
$\mathcal{C}=\mathcal{C}^{+}\cup\mathcal{C}^{-}$, where $\mathcal{C}^{+}$ runs
from $-\infty$ to $+\infty$, namely, forward in time, and $\mathcal{C}^{-}$
runs from $+\infty$ to $-\infty$, namely, backward in time. Physically, the
need of a backward and a forward in time direction is expected because, as
mentioned earlier, the density matrix is an operator so its time evolution
must include both time directions. On the Keldysh contour $\mathcal{C}$ we use
the definitions $\psi^{i}(t^{+})=\psi_{i}^{+}(t)$ and
$\psi^{i}(t^{-})=i\psi_{i}^{-}(t)$, where $\psi_{i}^{+}$ and $\psi_{i}^{-}$
are the independent Majoranas introduced in the action Eq. (10). Since the
action in the Keldysh contour is in general defined in the complex plane, it
is more natural to rewrite it as,

$\begin{split}iS=&-\int_{\mathcal{C}}dz\,\psi^{i}(z)\partial_{z}\psi^{i}(z)-i\int_{\mathcal{C}}dz\,i^{q/2}\sum_{i_{1}<\cdots<i_{q}}^{N}J_{i_{1}\cdots
i_{q}}\psi^{i_{1}}(z)\cdots\psi^{i_{q}}(z)\\\
&+\frac{\mu}{2}\int_{\mathcal{C}}dz\,dz^{\prime}\,K(z,z^{\prime})\psi^{i}(z)\psi^{i}(z^{\prime}),\end{split}$
(11)

where $dz=dt$ on $\mathcal{C}^{+}$, $dz=-dt$ on $\mathcal{C}^{-}$ and the
kernel $K(z,z^{\prime})$ is defined by

$\displaystyle\int_{\mathcal{C}}dz\,dz^{\prime}\,K(z,z^{\prime})\psi^{i}(z)\psi^{i}(z^{\prime})=\int
dt\left(\psi^{i}(t^{+})\psi^{i}(t^{-})-\psi^{i}(t^{-})\psi^{i}(t^{+})\right).$
(12)

For convenience, we did not include the irrelevant constant $-N\mu/2$. We note
that the resulting Keldysh contour (path) integral has no constraints related,
for instance, to initial conditions. This is another simplification that is
justified because we are not interested in the short time dynamics, namely,
implicitly we are assuming that the system will quickly forget about initial
conditions which is typical of quantum chaotic systems. We refer to Ref. [49]
for more details on the derivation of the action.

We now average over random couplings $J_{ijkl}$ and introduce two bi-local
functions $G^{ab}(z_{1},z_{2})$ and $\Sigma^{ab}(z_{1},z_{2})$ according to

$\displaystyle\delta(G^{ab}(t_{1},t_{2})+i\sum_{i}\psi^{a}_{i}(t_{1})\psi^{b}_{i}(t_{2}))=\int
d\Sigma
e^{-\frac{N}{2}\int\\!\\!\\!\int_{\mathcal{C}}\Sigma^{ab}(z_{1},z_{2})(G^{ab}(z_{1},z_{2})+i\sum_{i}\psi^{a}_{i}(z_{1})\psi^{b}_{i}(z_{2}))}$
(13)

with labels $a$ and $b$ $\in\\{+,-\\}$. After integrating the fermions, the
path integral $\int DGD\Sigma e^{iS}$ is a functional integral with respect to
$G$ and $\Sigma$ which can be evaluated by a saddle-point approximation. The
Schwinger-Dyson equations are given by (see Appendix C),

$\displaystyle(i\delta(z-z^{\prime})\partial_{z}-\Sigma(z,z^{\prime}))G(z,z^{\prime})$
$\displaystyle=$ $\displaystyle\mathrm{1}_{\mathcal{C}},$ (14)
$\displaystyle\Sigma(z,z^{\prime})+i^{q}J^{2}[G(z,z^{\prime})]^{q-1}-\frac{i}{2}\mu[K(z,z^{\prime})-K(z^{\prime},z)]$
$\displaystyle=$ $\displaystyle 0,$ (15)

where we have used that $G(z,z^{\prime})=-G(z^{\prime},z)$. In terms of the
$\pm$ components and the time variables, these equations can be written down
as

$\displaystyle\left(\begin{array}[]{cc}(\mathrm{i}\partial_{t}-\Sigma^{++})G^{++}+\Sigma^{+-}G^{-+}&(\mathrm{i}\partial_{t}-\Sigma^{++})G^{+-}+\Sigma^{+-}G^{--}\\\
-(\mathrm{i}\partial_{t}+\Sigma^{--})G^{-+}+\Sigma^{-+}G^{++}&-(\mathrm{i}\partial_{t}+\Sigma^{--})G^{--}+\Sigma^{-+}G^{+-}\end{array}\right)=\left(\begin{aligned}
\delta(t,t^{\prime})&\quad 0\\\
0\quad&\delta(t,t^{\prime})\end{aligned}\right),$ (18)

and

$\displaystyle\Sigma^{++}(t,t^{\prime})+J^{2}i^{q}G^{++}(t,t^{\prime})^{q-1}=0,$
$\displaystyle\Sigma^{--}(t,t^{\prime})+J^{2}i^{q}G^{--}(t,t^{\prime})^{q-1}=0,$
(19)
$\displaystyle\Sigma^{+-}(t,t^{\prime})+J^{2}i^{q}G^{+-}(t,t^{\prime})^{q-1}-i\mu\delta(t,t^{\prime})=0,$
$\displaystyle\Sigma^{-+}(t,t^{\prime})+J^{2}i^{q}G^{-+}(t,t^{\prime})^{q-1}+i\mu\delta(t,t^{\prime})=0.$

It is often convenient to express the solutions of the SD equations above in
terms of the real time Green’s functions

$\displaystyle G^{<}(t,t^{\prime})$ $\displaystyle=$ $\displaystyle\langle
G^{+-}(t,t^{\prime})\rangle,$ $\displaystyle G^{>}(t,t^{\prime})$
$\displaystyle=$ $\displaystyle\langle G^{-+}(t,t^{\prime})\rangle,$
$\displaystyle G^{T}(t,t^{\prime})$ $\displaystyle=$ $\displaystyle\langle
G^{++}(t,t^{\prime})\rangle,$ (20) $\displaystyle
G^{{\overline{T}}}(t,t^{\prime})$ $\displaystyle=$ $\displaystyle\langle
G^{--}(t,t^{\prime})\rangle.$ (21)

Similar combinations apply to the $\Sigma$ variables.

Because the Keldysh contour is time ordered, we have that

$\displaystyle G^{T}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\theta(t-t^{\prime})G^{>}(t,t^{\prime})+\theta(t^{\prime}-t)G^{<}(t,t^{\prime}),$
(22) $\displaystyle G^{\overline{T}}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\theta(t-t^{\prime})G^{<}(t,t^{\prime})+\theta(t^{\prime}-t)G^{>}(t,t^{\prime}).$
(23)

This gives the identity

$\displaystyle
G^{T}(t,t^{\prime})+G^{\overline{T}}(t,t^{\prime})=G^{>}(t,t^{\prime})+G^{<}(t,t^{\prime}).$
(24)

The action is invariant under $\psi(t^{+})\to\psi(t^{-})$,
$\psi(t^{-})\to\psi(t^{+})$ and $i\to-i$. Therefore,

$\displaystyle G^{>}(t,t^{\prime})=-G^{<}(t,t^{\prime}),$ (25)

or, in Fourier space, $G^{>}(\omega)=-G^{<}(\omega)$. Taking into account the
usual relation for finite temperature Green’s functions,
$G^{>}(\omega)=-\exp(\beta\omega)G^{<}(\omega)$, we conclude that the Green’s
functions are at infinite temperature, which is consistent with the infinite
temperature equilibrium state of the Lindblad evolution. In appendix B, we
will show that the conditions (24) and (25) also follow from the SD equations.

The Fourier transformed SD equations are given by

$\displaystyle(\omega-\Sigma^{T}(\omega))G^{T}(\omega)+\Sigma^{<}(\omega)G^{>}(\omega)$
$\displaystyle=$ $\displaystyle 1,$ (26)
$\displaystyle(\omega-\Sigma^{T}(\omega))G^{<}(\omega)+\Sigma^{<}(\omega)G^{\overline{T}}(\omega)$
$\displaystyle=$ $\displaystyle 0,$ (27)
$\displaystyle(\omega+\Sigma^{\overline{T}}(\omega))G^{>}(\omega)-\Sigma^{>}(\omega)G^{{T}}(\omega)$
$\displaystyle=$ $\displaystyle 0,$ (28)
$\displaystyle-(\omega+\Sigma^{\overline{T}}(\omega))G^{\overline{T}}(\omega)+\Sigma^{>}(\omega)G^{{<}}(\omega)$
$\displaystyle=$ $\displaystyle 1.$ (29)

Generally, the solution of the SD equations requires an $i\epsilon$
prescription, $\omega\to\omega+i\epsilon$, but for given boundary conditions,
a nonzero value of $\mu$ already selects solutions with the correct causal
structure.

In order to obtain the Lyapunov exponent, our first task is the evaluation of
the retarded Green function $G^{R}(t)=\theta(t)(G^{>}(t)-G^{<}(t))$. By using
that the solutions of the SD equations satisfy the conditions (24) and (25),
the third and fourth SD equations coincide with the first two. The difference
and the sum of the first two SD equations results in

$\displaystyle(\omega-\Sigma^{R}(\omega))G^{R}(\omega)=1,\qquad(\omega-\Sigma^{A}(\omega))G^{A}(\omega)=1.$
(30)

where we have defined the retarded and advanced Green’s function,

$\displaystyle G^{R}(\omega)$ $\displaystyle=$ $\displaystyle
G^{T}(\omega)-G^{<}(\omega),$ (31) $\displaystyle G^{A}(\omega)$
$\displaystyle=$ $\displaystyle G^{T}(\omega)+G^{<}(\omega),$ (32)

with analogous expressions for $\Sigma^{R,A}$. To solve the SD equations, we
will employ the reality conditions of the Green’s functions which we will
derive next. Following arguments in Ref. [83] for a closely related model,
they can also be obtained from the symmetry properties of the action.

The Euclidean time ordered Green’s function (denoted by a superscript $E$) in
the limit of interest is given by [78]

$\displaystyle
G^{E}(u,u^{\prime})=\left\langle\frac{1}{N}\Tr\sum_{i}\mathcal{T}\psi^{i}(u)\psi^{i}(u^{\prime})\right\rangle,$
(33)

with $\mathcal{T}$ the time ordering operator, and satisfies the analyticity
requirement [78]

$\displaystyle{G^{E}}^{*}(u,u^{\prime})=-G^{E}(-u^{*},-{u^{\prime}}^{*}).$
(34)

If the Hamiltonian is non-Hermitian, this identity is only valid if
$H^{\dagger}$ and $H$ have the same weight in the ensemble of random SYK
models. Therefore,

$\displaystyle{G^{>}}^{*}(\omega)$ $\displaystyle=$
$\displaystyle-{G^{>}}(\omega)={G^{<}}(\omega),$
$\displaystyle{G^{T}}^{*}(\omega)$ $\displaystyle=$ $\displaystyle
G^{T}(\omega)=-G^{{\overline{T}}}(\omega),$ (35)

and

$\displaystyle
G^{>}(\omega)=-G^{<}(\omega)=\frac{1}{2}({G^{R}(\omega)-G^{R}(\omega)^{*}}).$
(36)

We remind the reader that throughout this paper $\beta=0$. In terms of the
real time Green’s functions, the SD equations can be simplified to

$\displaystyle G^{<}(\omega)$ $\displaystyle=$
$\displaystyle-i{\mathrm{Im}G^{R}(\omega)},$
$\displaystyle\Sigma^{R}(\omega)+i\mu$ $\displaystyle=$
$\displaystyle-2iJ^{2}\int_{0}^{+\infty}dte^{i(\omega+i\epsilon)t}\mathrm{Im}[G^{<}(t)^{3}],$
$\displaystyle G^{R}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{(\omega+i\epsilon)-\Sigma^{R}(\omega)},$ (37)

where we have used that
$\Sigma^{R}(t)=\theta(t)(\Sigma^{>}(t)-\Sigma^{<}(t))$, Eq. (19),
$G^{>}(t)=(G^{<}(t))^{*}$, and $\int dt\delta(t)\theta(t)=1/2$. The real time
equations above are solved by standard iterative methods. Strictly speaking,
these expressions are only applicable to systems close to thermodynamical
equilibrium. As mentioned earlier, the steady state resulting from the
Lindbladian dynamics is [49] a TFD at infinite temperature $(\beta=0)$. At
sufficiently short times, the dynamics will depend on the initial state and
therefore cannot be described by the SD equations which were derived from the
Keldysh path integral without introducing any initial state constraint or
specifying a specific quench protocol. This independence on initial conditions
is anticipated, even for relatively short time scales, because classically
chaotic systems lose memory of the initial conditions rather quickly, so this
loss occurs much earlier than the exponential growth of the OTOC around the
Ehrenfest time due to the proliferation of quantum effects. However, for some
choices of initial conditions, the loss of memory may take longer to occur
which may lead to non-universal deviations from the OTOC results reported
below. In any case, we expect this approximation to be valid in most cases,
especially for initial states given by a TFD state at sufficiently high
temperature.

In the next section, we proceed with the computation of the OTOC around the
Ehrenfest time for an SYK coupled to a Markovian bath.

##  III Calculation of the OTOC

After solving the SD equations and finding an explicit expression for $G^{R}$,
we show in appendix C that the Lyapunov exponent $\lambda_{L}$ is computed by
solving

$\displaystyle\mathcal{F}(t_{1},t_{2})=\int
dt_{3}dt_{4}K_{R}(t_{1},t_{2};t_{3},t_{4})\mathcal{F}(t_{3},t_{4}),$ (38)

using the Ansatz

$\displaystyle\mathcal{F}(t_{1},t_{2})=e^{\lambda_{L}(t_{1}+t_{2})/2}f(t_{1}-t_{2}).$
(39)

The Bethe-Salpeter equation (38) is derived from the path integral formalism
by evaluating the leading $1/N$ correction and expressing it in terms of
Green’s functions resulting from the solution of the SD equations. It turns
out that the expression for the kernel $K_{R}(t_{1},t_{2};t_{3},t_{4})$ has
the same structure as the one for the single site SYK model [22],

$\displaystyle
K_{R}(t_{1},t_{2};t_{3},t_{4})=-(q-1)J^{2}G^{R}(t_{1}-t_{3})G^{R}(t_{2}-t_{4})G_{W}^{q-2}(t_{3}-t_{4}),$
(40)

where $G_{W}$ is the Wightman function, derived in appendix C, which can be
expressed as the analytic continuation of $G^{E}$, the Euclidean Green’s
function (33), $G_{W}(t)=G^{E}(it+\beta/2)$ with $\beta=0$ in this case.

It is in principle surprising that the structure of this kernel function is
the same as the one in the original SYK model [22] even though the SYK model
is now coupled to a Markovian bath. The reason is that the specific choice of
jump operators we are considering, linear in the Majorana fermions, leads to a
coupling term in the path integral which is linear in $G^{>}$ and $G^{<}$. As
a result, the leading perturbations to the large $N$ limit that are necessary
for the calculation of the OTOC only contribute linearly. Since the derivation
of the kernel function relies on second-order perturbations to the Green’s
functions, the form of the kernel function does not depend on the coupling
$\mu$, though the OTOC does depend on it through the retarded Green’s function
and the Wightman functions. More details can be found in appendix C. We stress
that this independence of the bath is a particularity of our choice of jump
operators. In general, jump operators involving two or more Majorana fermions
will also modify the structure of the kernel equation.

Using the Ansatz (39), the kernel equation (38) can be transformed into

$\displaystyle f(t_{1}-t_{2})$ $\displaystyle=$
$\displaystyle-(q-1)J^{2}\int_{-\infty}^{t_{1}}\int_{-\infty}^{t_{2}}dt_{3}dt_{4}e^{\frac{\lambda_{L}}{2}(t_{3}+t_{4}-t_{1}-t_{2})}$
$\displaystyle\times
G^{R}(t_{1}-t_{3})G^{R}(t_{2}-t_{4})G_{W}^{q-2}(t_{3}-t_{4})f(t_{3}-t_{4}).$

At time scales where the system is close to the equilibrium state, the
retarded Green’s function is expected to simplify to

$\displaystyle G_{R}(t)\approx-i\theta(t)e^{-\Gamma t},$ (42)

and the ladder equation reduces to a differential equation

$\displaystyle(\partial_{t_{1}}+\frac{\lambda_{L}}{2}+\Gamma)(\partial_{t_{2}}+\frac{\lambda_{L}}{2}+\Gamma)f(t_{1}-t_{2})=(q-1)J^{2}G_{W}^{q-2}(t_{1}-t_{2})f(t_{1}-t_{2}).$
(43)

Using translational invariance, this becomes the Schrödinger-like equation

$\displaystyle\left[-\partial_{t}^{2}+\left(\frac{\lambda_{L}}{2}+\Gamma\right)^{2}-(q-1)J^{2}G_{W}^{q-2}(t)\right]f(t)=0$
(44)

of a particle in the potential $-(q-1)J^{2}G_{W}^{q-2}(t)$. Note that
$G_{W}(t)$ is an even function of $t$. In the next section, we will
analytically solve this equation in the large $q$ limit. In the limit
$\beta\to 0$ of interest, we would have to solve the real time equations (37).
In that case, the Fourier transformed Wightman function is obtained from the
Euclidean Green’s function [78]

$\displaystyle G^{E}(u)=-\int\frac{d\omega}{2\pi}e^{-\omega
u}\mathrm{Im}[G^{R}(\omega)]$ (45)

resulting in the Wightman function $G^{W}(t)=G^{E}(it)$, and

$\displaystyle G_{W}(\omega)={\rm Im}[G^{R}(\omega)].$ (46)

The function $f(\omega)$ obeys the integral equation

$\displaystyle
f(\omega)=(q-1)J^{2}\left|G^{R}(\omega+i\frac{\lambda_{L}}{2})\right|^{2}\int\frac{d\omega^{\prime}}{2\pi}g_{W}(\omega-\omega^{\prime})f(\omega^{\prime}),$
(47)

where $g_{W}(\omega)$ is the Fourier transform of $G_{W}^{q-2}(t)$.
Numerically, after discretizing $\omega$, this equation becomes a finite
dimensional eigenvalue equation. Its eigenvalues depend on $\lambda_{L}$, and
$\lambda_{L}$ is determined by the condition that the integral equation has an
eigenvalue 1. Numerically, this can be done efficiently by bi-section in
$\lambda_{L}$.

Figure 1: The Lyapunov exponent $\lambda_{L}$ as a function of the coupling to
the bath $\mu$ from the numerical solution of the ladder equation (38) for
$q=4$. We observe that $\lambda_{L}$ decreases with $\mu$ and eventually
becomes negative for $\mu>0.33$.

###  III.1 Numerical OTOC and Lyapunov exponent

The SD equations are solved numerically by iterating the real-time equations
(37). After a Fourier transform, this gives the retarded Green’s function
$G^{R}(\omega)$. The Wightman function $G_{W}(\omega)$ is obtained from
$G^{R}(\omega)$ according to (46). This allows us to solve the kernel
equation, also called the Bethe-Salpeter equation (47), by discretizing the
frequencies and numerically solving the ensuing eigenvalue equation.

It was shown in Ref. [49] that for sufficiently long times and $\mu>0$, the
steady state resulting from the Liouvillian dynamics is a TFD state at
infinite temperature ($\beta=0$). Following the numerical procedure explained
above, we show in Fig. 1 the Lyapunov exponent $\lambda_{L}$ as a function of
the coupling to the bath for $q=4$. Interestingly, it shows a monotonous,
almost linear, decrease of the Lyapunov exponent with the coupling to the bath
$\mu$ until it becomes negative for $\mu>\mu_{c}\approx 0.33$. A negative
value of $\lambda_{L}$ indicates that the exponential growth of the OTOC,
which physically describes the increase in quantum uncertainty with time, has
been fully suppressed by the coupling to the bath. It is tempting to identify
this point as the transition between quantum chaotic dynamics, mostly driven
by the single-site SYK model, and non-chaotic dynamics mostly driven by the
environment. More work needs to be done to fully characterize the nature of
this transition.

###  III.2 Analytic calculation of $\lambda_{L}$ in the large $q$ limit

We now proceed to the analytic calculation of the Lyapunov exponent for
$\beta\to 0$ in the large $q$ limit at fixed

$\displaystyle\mathcal{J}^{2}=\frac{qJ^{2}}{2^{q-1}}\quad{\rm
and}\quad\widehat{\mu}=\mu q.$ (48)

In this limit, the self-energy becomes $O(1/q)$ which makes it possible to
solve the SD equations analytically. Following the derivation for the
Maldacena-Qi model [76], the real time Green’s functions can be expanded as,

$\displaystyle G^{T}(t)=-\frac{i}{2}\mathrm{{\rm
sgn}(t)}(1+\frac{1}{q}g^{T}(t)+\cdots),$ (49) $\displaystyle
G^{<}(t)=\frac{i}{2}(1+\frac{1}{q}g^{<}(t)+\cdots),$

and, for $\beta=0$, $G^{>}=-G^{<}$. Therefore,

$\displaystyle
G^{R}(t)=\theta(t)(G^{>}(t)-G^{<}(t))=-i\theta(t)(1+\frac{1}{q}g^{<}(t)+\cdots),$
(50)

and the self-energies are given by

$\displaystyle\Sigma^{T}(t)$ $\displaystyle=$
$\displaystyle-i^{q}J^{2}{G^{T}(t)}^{q-1}=\frac{i^{2q-1}J^{2}}{2^{q-1}}{\rm
sgn}^{q-1}(t)e^{g^{T}(t)},$ $\displaystyle\Sigma^{<}(t)$ $\displaystyle=$
$\displaystyle-i^{q}J^{2}{G^{<}(t)}^{q-1}+i\mu=\frac{i^{2q-1}J^{2}}{2^{q-1}}e^{g^{<}(t)}+i\mu.$
(51)

Substituting the large $q$ expressions for $G$ and $\Sigma$ into the
Schwinger-Dyson equations and collecting the $1/q$ contributions (the leading
order cancels) results in the differential equations

$\displaystyle\partial_{t}^{2}({\rm sgn}(t)g^{T}(t))+2\mathcal{J}^{2}{\rm
sgn}(t)e^{g^{T}(t)}=0,$
$\displaystyle\partial_{t}^{2}g^{<}(t)+2\mathcal{J}^{2}e^{g^{<}(t)}=0.$ (52)

The solutions of these equations are given by

$\displaystyle
e^{g^{T}(t)}=\frac{\alpha^{2}}{\mathcal{J}^{2}\cosh^{2}(\alpha|t|+\gamma)},\qquad
e^{g^{<}(t)}=\frac{\widetilde{\alpha}^{2}}{\mathcal{J}^{2}\cosh^{2}(\widetilde{\alpha}|t|+\widetilde{\gamma})},$
(53)

where $\alpha,\;\widetilde{\alpha},\;\gamma$ and $\widetilde{\gamma}$ are
integration constants determined by the conditions

$\displaystyle\frac{\alpha^{2}}{\mathcal{J}^{2}\cosh^{2}\gamma}=1,\quad
2\widetilde{\alpha}\tanh\widetilde{\gamma}=\widehat{\mu},\quad\widetilde{\alpha}=\alpha,\quad\widetilde{\gamma}=\gamma,$
(54)

resulting from the imposition of the following boundary conditions:
$G^{T}(0^{+})={i}/{2}$, so $\exp({g^{T}(0)})=1$,
$\partial_{t}g^{<}(0)=-\widehat{\mu}$, which can be derived by recalling the
Schwinger-Dyson equations and applying the ${1}/{q}$ expansion. The last two
conditions in (54) come from the relation $G^{T}(t)=-\mathrm{{\rm
sgn}}(t)G^{<}(t)$ in the $t\to\infty$ limit [76].

For large $t$, the retarded Green function $G^{R}(t)$ is expected to decay
exponentially, $G^{R}(t)\propto e^{-\Gamma t}$. This can be shown explicitly
in the limits of large $q$ and $t\gg 1$ [84, 76, 85, 49]. Using that the
Green’s function vanishes for $t\to\infty$ as

$\displaystyle G^{<}(t)=\frac{i}{2}e^{-2\alpha t/q},$ (55)

where the exponent is obtained by matching to the $t\ll q$ solution. Although
this result is only valid for $t\gg 1$, it also turns out to be a good
approximation for $t<1$ [86]. To enable a fully analytical calculation of the
Lyapunov exponent, we use this approximation so that the retarded Green’s
function is given by

$\displaystyle G^{R}(t)=-i\theta(t)e^{-2\alpha t/q}.$ (56)

The same result can be obtained by a na\̈mathrm{i}ve exponentiation of the
$1/q$ correction in (50). The accuracy of this approximation has been
confirmed by comparing with results from the numerical solution of the SD
equations at a large but finite $q$. We note that this approximation is
incompatible with a strict large $q$ expansion but, at finite $q$, as in the
Hermitian case [22], we shall see it gives a much better approximation. For
example, in the weak-coupling/high-temperature limit of the single site SYK
model, the $1/q$-correction to the Lyapunov exponent is given by
$2\mathcal{J}(1-2/q)$ so approximately for $q=4$ the Lyapunov exponent is
suppressed by a factor two in rough agreement with numerical results from
solving the SD equations (see Fig. 11 of [22]).

The $(q-2)$th power of the Wightman function (see Eq. (46) for an expression
of the Wightman function in terms of the retarded Green’s function) for
$\beta\to 0$ is given by,

$\displaystyle(q-1)J^{2}G_{W}^{q-2}(t)=(q-1)J^{2}(-i)^{q-2}{G^{<}}^{q-2}(t)=\frac{2\alpha^{2}}{\cosh^{2}(\alpha|t|+\gamma)}.$
(57)

Therefore, in the kernel equation (38), the main contributions to the integral
come from $t_{3}-t_{4}\sim O(1)$ and both $t_{1}-t_{3}$ and $t_{2}-t_{4}$
$\sim O(q)$.

From the discussion in the previous section, we know that for an exponentially
decreasing retarded Green’s function the kernel equation reduces to the
differential equation

$\displaystyle\left[-\partial_{t}^{2}+\left(\Gamma+\frac{\lambda_{L}}{2}\right)^{2}-\frac{2\alpha^{2}}{\cosh^{2}(\alpha|t|+\gamma)}\right]f(t)=0.$
(58)

This is a Schrödinger-like equation for which we are seeking bound state
solutions. Since the potential has a discontinuous derivative at $t=0$, we
have to separate solutions for $t<0$ (equivalent to $\gamma\to-\gamma$) and
$t>0$. In order to satisfy the differential equation, these two solutions must
be continuous and with a continuous derivative at $t=0$. The general solution
of this differential equation is given by,

$\displaystyle
f(t)=c_{1}e^{-\sigma(t+\gamma/\alpha)}(\sigma+\alpha\tanh(\alpha
t+\gamma))+c_{2}e^{\sigma(t-\gamma/\alpha)}(-\sigma+\alpha\tanh(\alpha
t-\gamma)).$ (59)

Using the boundary conditions that $f(|t|\to\infty)=0$, we have that $c_{2}=0$
for $t>0$ and $c_{1}=0$ for $t<0$ when also $\gamma\to-\gamma$. Continuity of
the solution at $t=0$ requires that $c_{2}=-c_{1}$. We thus find the solution

$\displaystyle f(t)$ $\displaystyle=$ $\displaystyle
c_{1}\theta(t)e^{-\sigma(t+\gamma/\alpha)}(\sigma+\alpha\tanh(\alpha(t+\gamma/\alpha))-c_{1}\theta(-t)e^{\sigma(t-\gamma/\alpha)}(-\sigma+\alpha\tanh(\alpha(t-\gamma/\alpha))$
(60) $\displaystyle=$ $\displaystyle
c_{1}e^{-\sigma(|t|)+\gamma/\alpha)}(\sigma+\alpha\tanh(\alpha|t|+\gamma)).$

The second derivative in Eq.(58) gives rise to a term
$f^{\prime}(t)\delta(t)$, which has to vanish for $f(t)$ to be a solution of
the differential equation (58), i.e. $f^{\prime}(0)=0$. Therefore,

$\displaystyle\frac{\sigma}{\alpha}\left(-\frac{\sigma}{\alpha}-\tanh\gamma)\right)+\frac{1}{\cosh^{2}\gamma}=0.$
(61)

Solving for ${\sigma}/{\alpha}$ leads to,

$\displaystyle\frac{\sigma}{\alpha}=-\frac{1}{2}\tanh\gamma+\frac{\sqrt{4+\sinh^{2}\gamma}}{2\cosh\gamma}.$
(62)

From Eq. (54) we obtain

$\displaystyle\sinh\gamma$ $\displaystyle=$
$\displaystyle\frac{\widehat{\mu}}{2\mathcal{J}},\qquad\alpha=\mathcal{J}\cosh\gamma.$
(63)

Combining these equations with the expression for $\sigma/\alpha$ (62), we can
solve for the Lyapunov exponent $\lambda_{L}$ and the decay rate $\Gamma$ as a
function of $\mu$,

$\displaystyle\lambda_{L}=2\sigma-\frac{4\alpha}{q},\qquad\Gamma\equiv\frac{2\alpha}{q}.$
(64)

Figure 2: Analytic Lyapunov exponent $\lambda_{L}$ as a function
$\widehat{\mu}=\mu q$ for various finite values of $q$ (black curves) and
$q\to\infty$ (red curve). It becomes negative for sufficiently large
${\widehat{\mu}}$.

.

In Fig. 2, we show the result of the Lyapunov exponent versus $\widehat{\mu}$
for various values of $q$. Beyond a critical value of $\mu$, the Lyapunov
exponent becomes negative which indicates the strong suppression of quantum
effects consistent with a transition from quantum chaotic to quasi-classical
dynamics. We note that although these are results from a large-$q$ analysis,
the Lyapunov exponent for $q=2$ still vanishes as was expected for a non
quantum chaotic system.

The resulting analytical expressions for $\Gamma$ and $\lambda_{L}$, together
with the numerical results for $q=16,\;24,\;32$ and $64$, from solving the SD
equations numerically, are shown in Fig. 3. The value of the coupling constant
is equal to $\mathcal{J}=\sqrt{q}{J}/{2^{{(q-1)}/{2}}}=1$. The agreement with
the analytical result improves for increasing values of $q$. Also when
$\lambda_{L}<0$, the SD equations and the Bethe-Salpeter equations can be
solved in the same way without any problem. The results show a similar
agreement with the analytical results, but we do not show this part of the
curves.

Figure 3: Lyapunov exponent $\lambda_{L}$ (left) and decay rate $\Gamma$
(right) for $\mathcal{J}=1$, $q=\leavevmode\nobreak\ 16,\leavevmode\nobreak\
32,\leavevmode\nobreak\ 48,\leavevmode\nobreak\ 64$, $\beta=0$ for various
$\widehat{\mu}$. The dots correspond to the numerical solution of the kernel
equation (47), while the solid curves stand for the analytic large $q$
prediction Eq. (64). We observe that the agreement becomes better as $q$
increases. They are already almost indistinguishable for $q=64$.

We have also checked that the solution, Eq. (60), satisfies the integral
equation (58) by an explicit evaluation of the integral (see appendix D). In
that case, the quantization condition (62) results from imposing that the
solution Eq. (60) is an eigenfunction of the integral kernel equation (58).

##  IV Discussion and Conclusions

The results of this paper shed light on the question of precisely defining
quantum chaos in dissipative systems. In principle, it is tempting to extend
the BGS conjecture to systems with a complex spectrum, namely to state that
dissipative quantum chaotic dynamics is characterized by the agreement of the
spectral correlations of the system, in our case an SYK model coupled to a
bath, with the predictions of an ensemble of non-Hermitian random matrices
with the corresponding symmetries. Indeed, at least in certain cases, it has
been found [39, 40] level statistics of non-Hermitian many-body quantum
systems are well described by non-Hermitian random matrix theory. However,
dissipation tends to suppress quantum effects, so it is rather unclear whether
the extension of the BGS conjecture to non-Hermitian systems is in all cases
meaningful, even if spectral correlations are well described by (non-
Hermitian) random matrix theory. The findings of this paper show that
dissipation suppresses the Lyapunov exponent and eventually leads to a
negative Lyapunov exponent for a sufficiently strong dissipation strength. In
that case, quantum fluctuations decrease with time, so the initial wave packet
stays close to the classical trajectory. This suggests that a necessary
condition for the existence of dissipative $quantum$ chaos is a positive
Lyapunov exponent. However, a warning is in order, we have studied a very
specific Markovian bath leading to a coupling term in the Liouvillian which is
both integrable and Hermitian. It would be necessary to confirm whether the
observed transition occurs for more general jump operators with a larger
number of Majorana fermions and random couplings [47, 84] leading to coupling
terms in the Liouvillian which are both non-Hermitian and potentially quantum
chaotic.

Another issue, also related to the definition of quantum chaos in Hermitian
systems, is the bound [13] on the Lyapunov exponent $\lambda_{L}\leq 2\pi
T/\hbar$ for systems at finite temperature $T$. The saturation of the bound is
a signature of many-body quantum chaotic systems with a gravity dual. It would
be interesting to explore the extension of the bound to dissipative quantum
chaos. Unfortunately, our setting leads to a steady state at infinite
temperature whose Lyapunov exponent is not restricted by the bound. Different
jump operators will generally lead to different steady states so fine tuning
may be necessary to equilibrate to a Gibbs state which complicates the study
of finite temperature effects.

Finally, we comment on the universality of our results. As was mentioned in
the introduction, a similar vanishing of the Lyapunov exponent has been
reported [75] in a two-site SYK model where one of the SYK models acts as a
thermal bath provided that the coupling to the system, the other SYK, is
strong enough. We do not believe that this result depends on the details of
the SYK model. Quite the contrary, to some extent, the SYK model is the most
quantum chaotic system due to the random couplings and the infinite range of
its interactions, so we would expect that interacting models with a sparser
interaction will be less robust to the suppression of quantum effects due to
the contact with a bath and the subsequent suppression of the exponential
growth of the OTOC when the coupling is sufficiently strong. However, as
mentioned earlier, additional results within the Lindbladian approach, using
more general jump operators, are certainly desirable to confirm this
expectation.

In conclusion, we have studied the OTOC for the vectorized formulation of a
single-site SYK model coupled to a Markovian bath described by the Lindblad
formalism. We have found that the growth is still exponential around the
Ehrenfest time, but the Lyapunov exponent decreases with the coupling strength
to the bath and eventually becomes negative leading to a transition whose
details remain to be studied. Therefore, a natural definition of quantum chaos
in dissipative systems is the existence of a positive Lyapunov exponent
independent of whether level statistics agrees with the corresponding non-
Hermitian random matrix prediction whose precise relation to dissipative
quantum chaos dynamics remains to be demonstrated. We note that a limitation
of this characterization is that it requires a well defined semiclassical
limit, in the SYK model it is the existence of a $1/N$ expansion. In certain
Hermitian systems, like certain spin-chains [87, 88, 89, 90] whose level
statistics is well described by random matrix theory [91], no exponential
growth of the OTOC is observed because of this reason.

As mentioned earlier, a natural extension of this work is to study the
dynamics for more general jump operators in order to clarify the details of
the transition. It would also be interesting to compute the Lyapunov exponent
in dissipative spin-chains and other similar settings in order to test the
universality of the suppression of chaos in dissipative systems. Moreover, it
would be worthwhile to extend our results to finite temperature and to work
out the gravity dual interpretation, especially the mentioned transition. We
note that our SYK setting has striking similarities with a wormhole
configuration in a global near de Sitter background in two dimensions [92, 93,
49, 94].

###### Acknowledgements.

We thank Lucas Sa for interesting discussions. A. M. G. G. acknowledges
illuminating correspondence with Victor Godet. A. M. G. G. and J. Z.
acknowledge support from the National Natural Science Foundation of China
(NSFC): Individual Grant No. 12374138, Research Fund for International Senior
Scientists No. 12350710180, and National Key R$\&$D Program of China (Project
ID: 2019YFA0308603). A. M. G. G. acknowledges support from a Shanghai talent
program. J. J. M. V. is supported in part by U.S. DOE Grant No. DE-
FAG-88FR40388.

##  Appendix A Conventions

In this appendix, we summarize some of the conventions we have been using in
this paper.

The Fourier transform and its inverse are defined by

$\displaystyle f(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega
t}f(\omega),$ (65) $\displaystyle f(\omega)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dte^{i\omega t}f(t).$ (66)

For translational invariant functions we write

$\displaystyle f(t_{1},t_{2})=f(t_{1}-t_{2})=f(t)$ (67)

.

##  Appendix B Solution of the Schwinger-Dyson Equations

In this appendix, we show that the solutions of the SD equations satisfy the
properites

$\displaystyle G^{>}(\omega)$ $\displaystyle=$ $\displaystyle-G^{<}(\omega),$
$\displaystyle G^{T}(\omega)$ $\displaystyle=$
$\displaystyle-G^{\overline{T}}(\omega).$ (68)

To solve the SD equations we introduce

$\displaystyle\widetilde{\Sigma}^{+-}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\Sigma^{+-}(t,t^{\prime})-i\mu\delta(t,t^{\prime}),$
$\displaystyle\widetilde{\Sigma}^{-+}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\Sigma^{-+}(t,t^{\prime})+i\mu\delta(t,t^{\prime}).$ (69)

Then the $\widetilde{\Sigma}^{\alpha\beta}$, with $\alpha,\beta\in\\{+,-\\}$,
have the same transformation properties under $t\to-t$ as the corresponding
Green’s functions. For the Fourier transformed Green’s functions and self-
energies, this results in

$\displaystyle\widetilde{\Sigma}^{>}(-\omega)=-\widetilde{\Sigma}^{<}(\omega).$
(70)

We thus find the SD equations

$\displaystyle(\omega-\Sigma^{T}(\omega))G^{T}(\omega)+(\widetilde{\Sigma}^{<}(\omega)+i\mu)G^{>}(\omega)$
$\displaystyle=$ $\displaystyle 1,$ (71)
$\displaystyle(\omega-\Sigma^{T}(\omega))G^{<}(\omega)+(\widetilde{\Sigma}^{<}(\omega)+i\mu)G^{\overline{T}}(\omega)$
$\displaystyle=$ $\displaystyle 0,$ (72)
$\displaystyle(\omega+\Sigma^{\overline{T}}(\omega))G^{>}(\omega)-(\widetilde{\Sigma}^{>}(\omega)-i\mu)G^{{T}}(\omega)$
$\displaystyle=$ $\displaystyle 0,$ (73)
$\displaystyle-(\omega+\Sigma^{\overline{T}}(\omega))G^{\overline{T}}(\omega)+(\widetilde{\Sigma}^{>}(\omega)-i\mu)G^{{<}}(\omega)$
$\displaystyle=$ $\displaystyle 1.$ (74)

Adding the second and third equation and subtracting the fourth equation from
the first one gives

$\displaystyle(\omega-\Sigma^{T}(\omega))G^{<}(\omega)+\widetilde{\Sigma}^{<}(\omega)G^{\overline{T}}(\omega)+(\omega+\Sigma^{\overline{T}}(\omega))G^{>}(\omega)-\widetilde{\Sigma}^{>}(\omega)G^{{T}}(\omega)$
$\displaystyle=$
$\displaystyle-i\mu(G^{\overline{T}}(\omega)+G^{{{\overline{T}}}}(\omega)),$
$\displaystyle(\omega-\Sigma^{T}(\omega))G^{T}(\omega)+\widetilde{\Sigma}^{<}(\omega)G^{>}(\omega)+(\omega+\Sigma^{\overline{T}}(\omega))G^{\overline{T}}(\omega)-\widetilde{\Sigma}^{>}(\omega)G^{{<}}(\omega)$
$\displaystyle=$ $\displaystyle i\mu(G^{>}(\omega)+G^{<}(\omega)).$

The combination of the previous equations results in

$\displaystyle\widetilde{\Sigma}^{<}(\omega)G^{>}(\omega)-\widetilde{\Sigma}^{>}(\omega)G^{{<}}(\omega)$
$\displaystyle=$ $\displaystyle i\mu(G^{>}(\omega)+G^{<}(\omega)),$ (77)

which can be rewritten as

$\displaystyle\frac{1}{2}(\widetilde{\Sigma}^{<}(\omega)-\widetilde{\Sigma}^{>}(\omega))(G^{<}(\omega)+G^{>}(\omega))-\frac{1}{2}(\widetilde{\Sigma}^{<}(\omega)+\widetilde{\Sigma}^{>}(\omega))(G^{<}(\omega)-G^{>}(\omega))$
$\displaystyle=$ $\displaystyle i\mu(G^{>}(\omega)+G^{<}(\omega)).$

After substituting the equations for the self-energy, we obtain one equation
with two unknown functions, resulting in a relation between $G^{>}(\omega)$
and $G^{<}(\omega)$. We expect at most a finite set of solutions. If we make
the Ansatz $G^{>}(\omega)=\alpha\exp(b\omega)G^{<}(\omega)$ we find

$\displaystyle\alpha
e^{b\omega}(1-\alpha^{2})\Sigma^{<}(\omega)G^{<}(\omega)=i\mu(1+\alpha
e^{b\omega})G^{<}(\omega).$ (79)

Since both $G^{<}(\omega)$ and $\Sigma^{<}(\omega)$ are non-vanishing, not
constant functions, we conclude that $\alpha=-1$ and $b=0$.

We thus have shown that the SD equations imply

$\displaystyle G^{>}(\omega)=-G^{<}(\omega).$ (80)

Using the usual finite temperature relation
$G^{>}(\omega)=-\exp(\beta\omega)G^{<}(\omega)$, we confirm the solutions of
the SD equations give Green’s functions at $\beta=0$ in agreement with the
equilibrium solution of the Lindblad equation equal to the identity, or in
vectorized form, the TFD state at $\beta=0$. After substituting (80) into (B)
we find

$\displaystyle-(\Sigma^{T}(\omega)+\Sigma^{T}(\omega))G^{<}(\omega)+\widetilde{\Sigma}^{<}(\omega)(G^{\overline{T}}(\omega)+G^{T}(\omega))=-i\mu(G^{\overline{T}}(\omega)+G^{{{\overline{T}}}}(\omega)).$
(81)

This shows that

$\displaystyle G^{\overline{T}}(\omega)=-G^{T}(\omega).$ (82)

##  Appendix C Derivation of the OTOC kernel equation

In this appendix, we derive an expression for the OTOC in terms of Green’s
functions. The derivation is along the lines of Ref. [22] and follows the
calculation of the OTOC [78] for the Maldacena-Qi model [76].

The regularized OTOC [22, 76] is defined by the four-point function

$\displaystyle
F(t_{1},t_{3})=\frac{1}{N^{2}}\frac{1}{Z}\left\langle\Tr\rho^{1/4}\psi_{i}^{+}(t_{1})\rho^{1/4}\psi_{j}^{-}(0)\rho^{1/4}\psi_{i}^{+}(t_{3})\rho^{1/4}\psi_{j}^{-}(0)\right\rangle,$
(83)

where $\psi_{i}^{\pm}(t)$ refers to Majoranas living on the forward and
backward in time sections of the contour Fig. 4 to be introduced shortly,
$\langle\cdots\rangle$ denote averaging over the random couplings and we have
introduced a “Boltzmann” factor $\rho=e^{-i\eta\mathcal{L}}$ as a regulator
for a proper definition of the Keldysh contour, see Fig. 4. We will take the
limit $\eta\to 0$ at the end of the calculation. The sum over repeated indexes
$i,j$ is assumed and the average is normalized by the ”partition” function
$Z=\Tr e^{-i\eta\mathcal{L}}$.

$u$$u$$u$$u$$\mathcal{C}$$0$$\frac{\eta}{4}$$\frac{\eta}{2}$$\frac{3\eta}{4}$$\eta$
Figure 4: In order to calculate the OTOC for the SYK model with dissipation,
we introduce the Keldysh contour $\mathcal{C}$ (solid line). The solid back
disks and the dashed lines denote the values of the arguments for which the
OTOC is evaluated.

In order to proceed, it is convenient to introduce the Euclidean four-point
function given by

$\displaystyle F^{E}(u_{1},u_{2},u_{3},u_{4})$
$\displaystyle=\frac{1}{N^{2}}\left\langle\Tr\psi_{i}(u_{1})\psi_{j}(u_{2})\psi_{i}(u_{3})\psi_{j}(u_{4})\right\rangle,$
(84)

where the expectation value is with respect to the path integral over Majorana
fermions discussed below and the Euclidean coordinates given by

$\displaystyle u_{1}=it_{1}+\frac{\eta}{4},\quad u_{2}=\frac{\eta}{2},\quad
u_{3}=it_{3}+\frac{3\eta}{4},\quad u_{4}=\eta,$ (85)

are time ordered on the Keldysh contour of Fig. 4. In this appendix we use the
Euclidean variable $u$ related to the contour integration variable $z$ in the
main text by $u=iz$. Below we will omit the redundant superscripts of the
$\psi_{i}$. The four-point function can be expressed as a two-replica path
integral that can be evaluated by a saddle-point approximation. To leading
order in $1/N$ the result factorizes into the product of two Green’s
functions, but these do not contribute to the exponentially increasing part of
the OTOC and we will ignore them. The exponentially increasing part of the
OTOC, ${\cal F}(t_{1},0,t_{3},0)$, is obtained from the connected $1/N$
corrections111There are also disconnected $1/N$ corrections from the
contraction of $\psi_{i}(u_{1})$ with $\psi_{j}(u_{4})$ and $\psi_{i}(u_{2})$
with $\psi_{j}(u_{3})$.

$\displaystyle
F(t_{1},0,t_{3},0)=-\frac{1}{N^{2}}\langle\psi_{i}(u_{1})\psi_{i}(u_{3})\rangle\langle\psi_{j}(u_{2})\psi_{j}(u_{4})\rangle+\frac{1}{N}{\cal
F}(u_{1},u_{2},u_{3},u_{4}),$ (86)

which is the average of the product of two fermion bilinears
$\sum_{k}\psi_{k}(u)\psi_{k}(v)$.

The four-point function is given by the expectation value Eq. (84) with
respect to the action

$\displaystyle-S$ $\displaystyle=$
$\displaystyle\int_{\mathcal{C}}\left(-\frac{i}{2}\sum\psi_{i}(u)\partial_{u}\psi_{i}(u)-i^{q/2+1}\sum
J_{ijkl}\psi_{i}(u)\psi_{j}(u)\psi_{k}(u)\psi_{l}(u)\right)$ (87)
$\displaystyle+\frac{\mu}{2}\int\\!\\!\\!\int_{\mathcal{C}}\sum_{i}\psi_{i}(u_{a})\psi_{i}(u_{b})K(u_{a},u_{b}),$

where $\mathcal{C}$ is the Keldysh contour shown in Fig. 4, and
$K(u,v)=-K(v,u)$ is defined by

$\displaystyle\frac{1}{2}\int\\!\\!\\!\int_{\mathcal{C}}\sum_{i}\psi_{i}(u_{a})\psi_{i}(u_{b})K(u_{a},u_{b})=\int_{0}^{\infty}dt\psi^{+}(t)\psi^{-}(t).$
(88)

The limit of integration of the integration contour $\mathcal{C}$ and the
right hand side of the above equation can be extended to $-\infty$ without
affecting the result. We note that to simplify the notation in this appendix,
we are not writing down the differentials $du$ for single integrals,
$du_{a}du_{b}$ for double integrals, and $du_{a}du_{b}du_{c}du_{d}$ for
quadruple integrals.

The corresponding partition function after averaging over the random couplings
is given by

$\displaystyle\langle Z\rangle_{J}$ $\displaystyle=\int
D\psi\exp\left(-\frac{i}{2}\int_{\mathcal{C}}\psi_{i}(u)\partial_{u}\psi_{i}(u)-\frac{N}{2q}J^{2}\int\\!\\!\\!\int_{\mathcal{C}}\left(\frac{1}{N}\psi_{i}(u_{a})\psi_{i}(u_{b})\right)^{q}+\frac{\mu}{2}\int_{\mathcal{C}}\psi_{i}(u_{a})\psi_{i}(u_{b})K(u_{a},u_{b})\right).$
(89)

Introducing the identity

$\displaystyle 1\sim\int
DGD\Sigma\exp\left(-\frac{N}{2}\int\\!\\!\\!\int_{\mathcal{C}}\Sigma(u_{a},u_{b})\left(G(u_{a},u_{b})+\frac{i}{N}\sum_{i}\psi_{i}(u_{a})\psi_{i}(u_{b})\right)\right)$
(90)

the partition function becomes (note that $q$ is even)

$\displaystyle\langle Z\rangle_{J}$ $\displaystyle=$ $\displaystyle\int
DGD\Sigma
D\psi\exp\left(-\frac{i}{2}\int_{\mathcal{C}}\psi_{i}(u)\partial_{u}\psi_{i}(u)-\frac{N}{2q}J^{2}\int\\!\\!\\!\int_{\mathcal{C}}\left[iG(u_{a},u_{b})\right]^{q}\right.$
$\displaystyle\left.-\frac{iN\mu}{2}\int\\!\\!\\!\int_{\mathcal{C}}G(u_{a},u_{b})K(u_{a},u_{b})-\frac{N}{2}\int\\!\\!\\!\int_{\mathcal{C}}\Sigma(u_{a},u_{b})G(u_{a},u_{b})-\frac{i}{2}\int\\!\\!\\!\int_{\mathcal{C}}\Sigma(u_{a},u_{b})\psi_{i}(u_{a})\psi_{i}(u_{b})\right).$

After integrating over the fermions, we obtain the $\Sigma G$ expression for
the partition function,

$\displaystyle\langle Z\rangle_{J}$ $\displaystyle=\int
DGD\Sigma\exp\left(\frac{N}{2}\int\int_{\mathcal{C}}\Tr\log\left(i\delta(u_{a},u_{b})\partial_{u_{a}}+i\Sigma(u_{a},u_{b})\right)-\frac{N}{2q}J^{2}\int\\!\\!\\!\int_{\mathcal{C}}\left[iG(u_{a},u_{b})\right]^{q}\right)$
(92)
$\displaystyle\times\exp\left(+\frac{iN\mu}{2}\int\\!\\!\\!\int_{\mathcal{C}}G(u_{a},u_{b})K(u_{a},u_{b})-\frac{N}{2}\int\\!\\!\\!\int_{\mathcal{C}}\Sigma(u_{a},u_{b})G(u_{a},u_{b})\right).$

Variation with respect to $\delta G(u_{a},u_{b})$ and
$\delta\Sigma(u_{a},u_{b})$ gives the Schwinger-Dyson equations

$\displaystyle-\Sigma(u_{a},u_{b})-i^{q}J^{2}G(u_{a},u_{b})^{q-1}+i\mu
K(u_{a},u_{b})\delta(u_{a},u_{b})=0,$ (93)
$\displaystyle(-\partial_{u}-\Sigma)\cdot G=1.$

In terms of the $z$ variables along the Keldysh contour of the main text, the
second equation becomes

$\displaystyle(i\partial_{z}-\Sigma)\cdot G=1,$ (94)

while the first equation remains the same. The OTOC is given by the leading
$1/N$ corrections, $\mathcal{F}(u_{1},u_{2},u_{3},u_{4})/N$, of

$\displaystyle F(t_{1},t_{2}=0,t_{3},t_{4}=0)=-\frac{1}{Z}\langle
G(u_{1},u_{3})G(u_{2},u_{4})\rangle$ (95)

with $u_{k}$ defined in Eq. (85). We will evaluate this correlation function
for general Euclidean $u_{k}$. The $1/N$ corrections are obtained from the
second order expansion of the action about the saddle point. Writing
$\Sigma(u,v)=\overline{\Sigma}(u,v)+\delta\Sigma(u,v)$ and
$G(u,v)=\overline{G}(u,v)+\delta G(u,v)$ we find

$\displaystyle-\frac{N}{4}\int\\!\\!\\!\int\\!\\!\\!\int\\!\\!\\!\int_{\mathcal{C}}\overline{G}(u_{a},u_{b})\delta\Sigma(u_{b},u_{c})\overline{G}(u_{c},u_{d})\delta\Sigma(u_{d},u_{a})-\frac{N}{2}\int\\!\\!\\!\int_{\mathcal{C}}\delta\Sigma(u_{a},u_{b})\delta
G(u_{a},u_{b})$
$\displaystyle-\frac{N}{4}J^{2}(q-1)i^{q}\int\\!\\!\\!\int_{\mathcal{C}}\overline{G}(u_{a},u_{b})^{q-2}\delta
G(u_{a},u_{b})^{2}.$ (96)

The integral over $\delta\Sigma(u_{a},u_{b})$ can be performed by completing
squares resulting in the effective action

$\displaystyle S^{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{N}{4}\int\\!\\!\\!\int\\!\\!\\!\int\\!\\!\\!\int_{C}\overline{G}^{-1}(u_{c},u_{d})\delta
G(u_{d},u_{a})\overline{G}^{-1}(u_{a},u_{b})\delta G(u_{b},u_{c})$ (97)
$\displaystyle-\frac{N}{4}J^{2}(q-1)i^{q}\int\\!\\!\\!\int_{\mathcal{C}}\overline{G}(u_{a},u_{b})^{q-2}\delta
G(u_{a},u_{b})^{2}$ $\displaystyle=$
$\displaystyle\frac{N}{4}\int\\!\\!\\!\int\\!\\!\\!\int\\!\\!\\!\int_{\mathcal{C}}\overline{G}^{-1}(u_{a},u_{c})\overline{G}^{-1}(u_{b},u_{d})\left[\delta(u_{a}-u_{f})\delta(u_{b}-u_{e})\right.$
$\displaystyle+\frac{N}{4}J^{2}(q-1)i^{q}\left.\overline{G}(u_{f},u_{a})\overline{G}(u_{e},u_{b})\overline{G}^{q-2}(u_{e},u_{f})\right]\delta
G(u_{e},u_{c})\delta G(u_{f},u_{d}).$

The $1/N$ correction to the four-point function

$\displaystyle{\cal F}(u_{1},u_{2},u_{3},u_{4})=-\int d[\delta G]\delta
G(u_{1},u_{3})\delta G(u_{2},u_{4})e^{-S^{\rm eff}}$ (98)

is evaluated as

$\displaystyle{\cal F}(u_{1},u_{2},u_{3},u_{4})$ $\displaystyle=$
$\displaystyle-\left[\mathbbm{1}-K\right]^{-1}(u_{1},u_{3};u,v)\overline{G}(u,u_{2})\overline{G}(v,u_{4})$
$\displaystyle=$
$\displaystyle-\left[\overline{G}(u_{3},u_{2})\overline{G}(u_{1},u_{4})-J^{2}(q-1)i^{q}\int\\!\\!\\!\int
K(u_{1},u_{3},u,v)\overline{G}(u,u_{2})\overline{G}(v,u_{4})+\cdots\right]$

with kernel $K$ given by

$\displaystyle
K(u_{a},u_{b},u_{c},u_{d})=-J^{2}(q-1)i^{q}\overline{G}(u_{a},u_{c})\overline{G}(u_{b},u_{d})\overline{G}^{q-2}(u_{c},u_{d}),$
(100)

and the points $u_{k}$ are equal to

$\displaystyle u_{1}=it_{1}+\frac{\eta}{4},\quad u_{2}=+\frac{\eta}{2},\quad
u_{3}=it_{3}+\frac{3\eta}{4},\quad u_{4}=+\eta.$ (101)

The action of the kernel $K$ by expanding the inverse in a geometric series
results in the combinations

$\displaystyle\int\\!\\!\\!\int_{\mathcal{C}}dudvK(it_{1}+\frac{\eta}{4},it_{3}+\frac{3\eta}{4},u,v)F(u,v).$
(102)

As discussed in Ref. [95, 78], the part of the $u$ integration that is on the
$3\eta/4$ branch cancels because the upward integration differs by a minus
sign from the downward integration (the $\epsilon$ difference between the
upward and downward path is negligible with respect to $3\eta/4$ because we
take the $\epsilon\to 0$ limit first). This is not the case when the $u$
integration is over the $\eta/4$ branch. Then, when $t_{1}$ is on the + path,

$\displaystyle
i(G(t_{1}^{+},t_{u}^{+})-G(t_{1}^{+},t_{u}^{-}))=i(G^{T}(t_{1},t_{u})-G^{<}(t_{1},t_{u}))=iG^{R}(t_{1},t_{u}).$
(103)

The result is the same when $t_{1}$ is on the minus path

$\displaystyle
i(G(t_{1}^{-},t_{u}^{+})-G(t_{1}^{-},t_{u}^{-}))=i(G^{>}(t_{1},t_{u})-G^{\overline{T}}(t_{1},t_{u}))=iG^{R}(t_{1},t_{u}).$
(104)

This also means that the $u$ variable only contributes when it is on the
$\eta/4$ branch. An analogue argument shows that the $v$ integration only
contributes when it is on the $3\eta/4$ branch. Combining the $+$ and $-$
paths again gives the retarded Green’s function $G^{R}(t_{3},t_{v}).$ In the
domain where $u$ and $v$ give non-vanishing contributions, they are separated
by $\eta/2$. Therefore the Green’s function $G(u,v)$ simply becomes the
Wightman function $G(it_{u},it_{v}+\eta/2)\equiv
G_{W}(i(t_{u}-t_{b})-\eta/2)$.

We conclude that the kernel is given by

$\displaystyle
K(t_{1},t_{2},t_{3},t_{4})=-J^{2}(q-1)G^{R}(t_{1},t_{3})G^{R}(t_{2},t_{4})G_{W}^{q-2}(t_{3},t_{4}).$
(105)

In the domain of exponential growth, for large $t_{1}$ and $t_{2}$, we have
that $(\mathbbm{1}-K)\mathcal{F}\ll\mathcal{F}_{0}$ with
$\mathcal{F}_{0}(u_{1},u_{2},u_{3},u_{4})=\overline{G}(u_{1},u_{3})\overline{G}(u_{2},u_{4})$.
Neglecting the right hand side, the Lyapunov exponent is obtained by solving
the integral equation

$\displaystyle(\mathbbm{1}-K)\mathcal{F}=0.$ (106)

##  Appendix D Solution of kernel equation for $\beta=0$

In this Appendix, we calculate the integral in Eq. (38) of the main text,

$\displaystyle I=\iint
dt_{3}dt_{4}\theta(t_{1}-t_{3})\theta(t_{2}-t_{4})e^{-\sigma(t_{1}-t_{3})}e^{-\sigma(t_{2}-t_{4})}\frac{2\alpha^{2}}{\cosh^{2}(\alpha|t_{3}-t_{4}|+\gamma)}f(|t_{3}-t_{4})|),$
(107)

where $f(|t_{3}-t_{4})|)$ is the solution of the differential equation (58),

$\displaystyle(\partial_{t}^{2}-\sigma^{2})f(t)=-\frac{2\alpha^{2}}{\cosh^{2}(\alpha|t|+\gamma)}f(t).$
(108)

We will check that the integral (107) is given by $f(t_{12})$. The integral
can be carried out by changing coordinates as follows,

$\displaystyle s=t_{3}-t_{4}\quad{\rm and}\quad t=\frac{1}{2}(t_{3}+t_{4}).$
(109)

This results in

$\displaystyle I=\int_{0}^{\infty}dte^{-2\sigma
t}\int_{-2t}^{2t}ds\frac{2\alpha^{2}}{\cosh^{2}(\alpha|s+t_{12}|+\gamma)}f(|s+t_{12}|),$
(110)

where $t_{12}=t_{1}-t_{2}$. After a partial integration of the $t$ integral,
integrating over $s$ results in

$\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{1}{2\sigma}\int_{0}^{\infty}dte^{-2\sigma
t}\left(\frac{2\alpha^{2}}{\cosh^{2}(\alpha|2t+t_{12}|+\gamma)}f(|2t+t_{12}|)+\frac{2\alpha^{2}}{\cosh^{2}(\alpha|-2t+t_{12}|+\gamma)}f(|-2t+t_{12}|)\right)$
(111) $\displaystyle=$
$\displaystyle\frac{1}{\sigma}\int_{t_{12}}^{\infty}due^{-\sigma(u-t_{12})}\frac{2\alpha^{2}}{\cosh^{2}(\alpha|u|+\gamma)}f(|u|)+\frac{1}{\sigma}\int^{t_{12}}_{-\infty}due^{\sigma(u-t_{12})}\frac{2\alpha^{2}}{\cosh^{2}(\alpha|u|+\gamma)}f(|u|).$

Next we use that $f(u)$ satisfies the differential equation (108) resulting in

$\displaystyle I$ $\displaystyle=$
$\displaystyle\frac{1}{2\sigma}\int_{t_{12}}^{\infty}due^{-\sigma(u-t_{12})}(\sigma^{2}-\partial_{u}^{2})f(|u|)+\frac{1}{2\sigma}\int^{t_{12}}_{0}due^{\sigma(u-t_{12})}(\sigma^{2}-\partial_{u}^{2})f(|u|)$
(112)
$\displaystyle+\frac{1}{2\sigma}\int^{0}_{-\infty}due^{\sigma(u-t_{12})}(\sigma^{2}-\partial_{u}^{2})f(|u|).$

The $\sigma^{2}$ term cancels against the integral that remains after two
partial integrations of the $\partial_{u}^{2}$ term leaving us with only
boundary contributions given by

$\displaystyle I$ $\displaystyle=$
$\displaystyle-\frac{1}{2\sigma}\left.e^{-\sigma(u-t_{12})}\partial_{u}f(|u|)\right|_{t_{12}}^{\infty}-\frac{1}{2}\left.e^{-\sigma(u-t_{12})}f(|u|)\right|_{t_{12}}^{\infty}$
(113)
$\displaystyle-\frac{1}{2\sigma}\left.e^{\sigma(u-t_{12})}\partial_{u}f(|u|)\right|^{t_{12}}_{0_{+}}-\frac{1}{2}\left.e^{\sigma(u-t_{12})}f(|u|)\right|^{t_{12}}_{0_{+}}$
$\displaystyle-\frac{1}{2\sigma}\left.e^{-\sigma
t_{12}}\partial_{u}f(|u|)\right|_{-\infty}^{0_{-}}+\frac{1}{2}\left.e^{-\sigma
t_{12}}f(|u|)\right|_{-\infty}^{0_{-}}.$

The solution $f(|u|)$ is a function of $|u|$ so that the derivatives at
$0_{-}$ and at $0_{+}$ have the opposite sign. We finally obtain

$\displaystyle
I=f(t_{12})+\frac{1}{2\sigma}\left.e^{\sigma(-t_{12})}\partial_{u}f(|u|)\right|_{0_{-}}^{0_{+}}.$
(114)

Therefore, the integral $I=f(t_{12})$ if $f^{\prime}(u=0)=0$, which is the
same condition as the one obtained from solving the differential equation in
the main text.

## References

  * Sekino and Susskind [2008] Y. Sekino and L. Susskind, Fast scramblers, Journal of High Energy Physics 10, 065 (2008).
  * Bagrets  _et al._ [2017] D. Bagrets, A. Altland, and A. Kamenev, Power-law out of time order correlation functions in the SYK model, Nuclear Physics B 921, 727 (2017).
  * Larkin and Ovchinnikov [1969] A. Larkin and Y. N. Ovchinnikov, Quasiclassical method in the theory of superconductivity, Sov Phys JETP 28, 1200 (1969).
  * Berman and Zaslavsky [1978] G. Berman and G. Zaslavsky, Condition of stochasticity in quantum nonlinear systems, Physica A: Statistical Mechanics and its Applications 91, 450 (1978).
  * Bohigas  _et al._ [1984] O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of Chaotic Quantum Spectra and Universality of Level Fluctuation Laws, Phys. Rev. Lett. 52, 1 (1984).
  * Wigner [1951] E. Wigner, On the statistical distribution of the widths and spacings of nuclear resonance levels, Math. Proc. Cam. Phil. Soc. 49, 790 (1951).
  * Dyson [1962a] F. J. Dyson, Statistical theory of the energy levels of complex systems. I, J. Math. Phys. 3, 140 (1962a).
  * Dyson [1962b] F. Dyson, Statistical Theory of the Energy Levels of Complex Systems. II, J. Math. Phys. 3, 157 (1962b).
  * Dyson [1962c] F. Dyson, Statistical Theory of the Energy Levels of Complex Systems. III, J. Math. Phys. 3, 166 (1962c).
  * Dyson [1962d] F. Dyson, The Threefold Way. Algebraic Structure of Symmetry Groups and Ensembles in Quantum Mechanics, J. Math. Phys. 3, 1199 (1962d).
  * Dyson [1972] F. Dyson, A Class of Matrix Ensembles, J. Math. Phys. 13, 90 (1972).
  * Mehta [2004] M. L. Mehta,  _Random matrices_ (Academic press, 2004).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 08, 106 (2016), arXiv:1503.01409 [hep-th] .
  * Kitaev [2015] A. Kitaev, A simple model of quantum holography (2015), string seminar at KITP and Entanglement 2015 program, 12 February, 7 April and 27 May 2015, http://online.kitp.ucsb.edu/online/entangled15/.
  * French and Wong [1970] J. French and S. Wong, Validity of random matrix theories for many-particle systems, Physics Letters B 33, 449 (1970).
  * Bohigas and Flores [1971a] O. Bohigas and J. Flores, Two-body random hamiltonian and level density, Physics Letters B 34, 261 (1971a).
  * French and Wong [1971] J. French and S. Wong, Some random-matrix level and spacing distributions for fixed-particle-rank interactions, Physics Letters B 35, 5 (1971).
  * Bohigas and Flores [1971b] O. Bohigas and J. Flores, Spacing and individual eigenvalue distributions of two-body random hamiltonians, Physics Letters B 35, 383 (1971b).
  * Benet  _et al._ [2001] L. Benet, T. Rupp, and H. A. Weidenmüller, Nonuniversal behavior of the $\mathit{k}$-body embedded gaussian unitary ensemble of random matrices, Phys. Rev. Lett. 87, 010601 (2001), arXiv:cond-mat/0010425 [cond-mat] .
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Gapless spin-fluid ground state in a random quantum Heisenberg magnet, Phys. Rev. Lett. 70, 3339 (1993), arXiv:cond-mat/9212030 [cond-mat] .
  * Sachdev [2010] S. Sachdev, Holographic Metals and the Fractionalized Fermi Liquid, Phys. Rev. Lett. 105, 151602 (2010).
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the Sachdev-Ye-Kitaev model, Phys. Rev. D 94, 106002 (2016).
  * Kobrin  _et al._ [2021] B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, Many-Body Chaos in the Sachdev-Ye-Kitaev Model, Phys. Rev. Lett. 126, 030602 (2021), arXiv:2002.05725 [hep-th] .
  * Garc\́mathrm{missing}ia-Garc\́mathrm{missing}ia  _et al._ [2023] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, C. Liu, and J. J. M. Verbaarschot, Sparsity independent Lyapunov exponent in the Sachdev-Ye-Kitaev model (2023), arXiv:2311.00639 [hep-th] .
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a and Verbaarschot [2016] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a and J. J. M. Verbaarschot, Spectral and thermodynamic properties of the Sachdev-Ye-Kitaev model, Phys. Rev. D 94, 126010 (2016), arXiv:1610.03816 [hep-th] .
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real spectra in non-hermitian hamiltonians having $PT$ symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * Wiersig [2014] J. Wiersig, Enhancing the sensitivity of frequency and energy splitting detection by using exceptional points: Application to microcavity sensors for single-particle detection, Phys. Rev. Lett. 112, 203901 (2014).
  * Gu  _et al._ [2020] Y. Gu, A. Kitaev, S. Sachdev, and G. Tarnopolsky, Notes on the complex Sachdev-Ye-Kitaev model, Journal of High Energy Physics 2020, 10.1007/jhep02(2020)157 (2020).
  * Kawabata  _et al._ [2018] K. Kawabata, Y. Ashida, H. Katsura, and M. Ueda, Parity-time-symmetric topological superconductor, Phys. Rev. B 98, 085116 (2018).
  * Kawabata  _et al._ [2017] K. Kawabata, Y. Ashida, and M. Ueda, Information retrieval and criticality in parity-time-symmetric systems, Phys. Rev. Lett. 119, 190401 (2017).
  * Li  _et al._ [2019] J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Observation of parity-time symmetry breaking transitions in a dissipative floquet system of ultracold atoms, Nature Communications 10, 10.1038/s41467-019-08596-1 (2019).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a and Godet [2021] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a and V. Godet, Euclidean wormhole in the Sachdev-Ye-Kitaev model, Phys. Rev. D 103, 046014 (2021), arXiv:2010.11633 [hep-th] .
  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Non-hermitian physics, Advances in Physics 69, 249 (2020).
  * Grobe  _et al._ [1988] R. Grobe, F. Haake, and H.-J. Sommers, Quantum Distinction of Regular and Chaotic Dissipative Motion, Phys. Rev. Lett. 61, 1899 (1988).
  * Fyodorov  _et al._ [1997] Y. V. Fyodorov, B. A. Khoruzhenko, and H.-J. Sommers, Almost Hermitian Random Matrices: Crossover from Wigner-Dyson to Ginibre Eigenvalue Statistics, Physical Review Letters 79, 557 (1997), cond-mat/9703152 .
  * Sá  _et al._ [2020] L. Sá, P. Ribeiro, and T. Prosen, Spectral and steady-state properties of random Liouvillians, J. Phys. A: Math. Theor. 53, 305303 (2020).
  * Sá  _et al._ [2020] L. Sá, P. Ribeiro, T. Can, and T. Prosen, Spectral transitions and universal steady states in random Kraus maps and circuits, Phys. Rev. B 102, 134310 (2020).
  * Rubio-Garc\́mathrm{i}a  _et al._ [2022] A. Rubio-Garc\́mathrm{i}a, R. Molina, and J. Dukelsky, From integrability to chaos in quantum liouvillians, SciPost Physics Core 5, 10.21468/scipostphyscore.5.2.026 (2022).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a  _et al._ [2022a] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, L. Sá, and J. J. M. Verbaarschot, Symmetry Classification and Universality in Non-Hermitian Many-Body Quantum Chaos by the Sachdev-Ye-Kitaev Model, Phys. Rev. X 12, 021040 (2022a).
  * Li  _et al._ [2021] J. Li, T. Prosen, and A. Chan, Spectral statistics of non-hermitian matrices and dissipative quantum chaos, Phys. Rev. Lett. 127, 170602 (2021).
  * Altland  _et al._ [2021] A. Altland, M. Fleischhauer, and S. Diehl, Symmetry classes of open fermionic quantum matter, Physical Review X 11, 021037 (2021).
  * Costa  _et al._ [2023] J. Costa, P. Ribeiro, A. De Luca, T. Prosen, and L. Sá, Spectral and steady-state properties of fermionic random quadratic liouvillians, SciPost Physics 15, 10.21468/scipostphys.15.4.145 (2023).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-hermitian physics, Phys. Rev. X 9, 041015 (2019), arXiv:1812.09133 [cond-mat.mes-hall] .
  * Ginibre [1965] J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6, 440 (1965).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a  _et al._ [2002] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, S. M. Nishigaki, and J. J. M. Verbaarschot, Critical statistics for non-hermitian matrices, Phys. Rev. E 66, 016132 (2002).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a  _et al._ [2023a] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, L. Sá, and J. J. M. Verbaarschot, Universality and its limits in non-hermitian many-body quantum chaos using the Sachdev-Ye-Kitaev model, Phys. Rev. D 107, 066007 (2023a).
  * Sá  _et al._ [2022] L. Sá, P. Ribeiro, and T. Prosen, Lindbladian dissipation of strongly-correlated quantum matter, Phys. Rev. Research 4, L022068 (2022).
  * Kulkarni  _et al._ [2022] A. Kulkarni, T. Numasawa, and S. Ryu, Lindbladian dynamics of the Sachdev-Ye-Kitaev model, Phys. Rev. B 106, 075138 (2022).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a  _et al._ [2023b] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, L. Sá, J. J. M. Verbaarschot, and J. P. Zheng, Keldysh wormholes and anomalous relaxation in the dissipative Sachdev-Ye-Kitaev model, Phys. Rev. D 107, 106006 (2023b).
  * Belavin  _et al._ [1969] A. A. Belavin, B. Y. Zeldovich, A. M. Perelomov, and V. S. Popov, Relaxation of quantum systems with equidistant spectra, Sov. Phys. JETP 29, 145 (1969).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of $N$-level systems, J. Math. Phys. 17, 821 (1976).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Manzano [2020] D. Manzano, A short introduction to the Lindblad master equation, AIP Advances 10, 025106 (2020).
  * Bergamasco  _et al._ [2023] P. D. Bergamasco, G. G. Carlo, and A. M. Rivas, Quantum lyapunov exponent in dissipative systems, Physical Review E 108, 024208 (2023).
  * Yoshida and Yao [2019] B. Yoshida and N. Y. Yao, Disentangling scrambling and decoherence via quantum teleportation, Physical Review X 9, 011006 (2019).
  * Tuziemski [2019] J. Tuziemski, Out-of-time-ordered correlation functions in open systems: A Feynman-Vernon influence functional approach, Physical Review A 100, 062106 (2019).
  * Weinstein  _et al._ [2023] Z. Weinstein, S. P. Kelly, J. Marino, and E. Altman, Scrambling transition in a radiative random unitary circuit, Physical Review Letters 131, 220404 (2023).
  * Zanardi and Anand [2021] P. Zanardi and N. Anand, Information scrambling and chaos in open quantum systems, Physical Review A 103, 062214 (2021).
  * qian Huang  _et al._ [2020] K. qian Huang, J. Wang, W.-L. Zhao, and J. Liu, Chaotic dynamics of a non-hermitian kicked particle, Journal of Physics: Condensed Matter 33, 055402 (2020).
  * Syzranov  _et al._ [2018] S. Syzranov, A. Gorshkov, and V. Galitski, Out-of-time-order correlators in finite open systems, Physical Review B 97, 161114 (2018).
  * Zhai and Yin [2020] L.-J. Zhai and S. Yin, Out-of-time-ordered correlator in non-hermitian quantum systems, Phys. Rev. B 102, 054303 (2020).
  * Turiaci [2019] G. J. Turiaci, An inelastic bound on chaos, Journal of High Energy Physics 2019, 10.1007/jhep07(2019)099 (2019).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a  _et al._ [2022b] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, V. Godet, C. Yin, and J. P. Zheng, Euclidean-to-Lorentzian wormhole transition and gravitational symmetry breaking in the Sachdev-Ye-Kitaev model, Phys. Rev. D 106, 046008 (2022b).
  * Schuster and Yao [2023] T. Schuster and N. Y. Yao, Operator Growth in Open Quantum Systems, Phys. Rev. Lett. 131, 160402 (2023), arXiv:2208.12272 [quant-ph] .
  * Bhattacharya  _et al._ [2022] A. Bhattacharya, P. Nandy, P. P. Nath, and H. Sahu, Operator growth and Krylov construction in dissipative open quantum systems, Journal of High Energy Physics 2022, 10.1007/jhep12(2022)081 (2022).
  * Bhattacharjee  _et al._ [2023a] B. Bhattacharjee, X. Cao, P. Nandy, and T. Pathak, Operator growth in open quantum systems: lessons from the dissipative SYK, Journal of High Energy Physics 2023, 10.1007/jhep03(2023)054 (2023a).
  * Zhang  _et al._ [2021] P. Zhang, S.-K. Jian, C. Liu, and X. Chen, Emergent Replica Conformal Symmetry in Non-Hermitian SYK2 Chains, Quantum 5, 579 (2021).
  * Jian  _et al._ [2021] S.-K. Jian, C. Liu, X. Chen, B. Swingle, and P. Zhang, Measurement-Induced Phase Transition in the Monitored Sachdev-Ye-Kitaev Model, Phys. Rev. Lett. 127, 140601 (2021), arXiv:2106.09635 [quant-ph] .
  * Bhattacharjee  _et al._ [2023b] B. Bhattacharjee, P. Nandy, and T. Pathak, Operator dynamics in Lindbladian SYK: a Krylov complexity perspective (2023b), arXiv:2311.00753 [quant-ph] .
  * Liu  _et al._ [2023] C. Liu, H. Tang, and H. Zhai, Krylov complexity in open quantum systems, Phys. Rev. Res. 5, 033085 (2023).
  * Srivatsa and von Keyserlingk [2023] N. S. Srivatsa and C. von Keyserlingk, The operator growth hypothesis in open quantum systems (2023), arXiv:2310.15376 [quant-ph] .
  * Zhang and Yu [2023] P. Zhang and Z. Yu, Dynamical transition of operator size growth in quantum systems embedded in an environment, Phys. Rev. Lett. 130, 250401 (2023).
  * Liu  _et al._ [2024] J. Liu, R. Meyer, and Z.-Y. Xian, Operator size growth in Lindbladian SYK (2024), arXiv:2403.07115 [hep-th] .
  * Chen  _et al._ [2017] Y. Chen, H. Zhai, and P. Zhang, Tunable quantum chaos in the Sachdev-Ye-Kitaev model coupled to a thermal bath, Journal of High Energy Physics 2017, 10.1007/jhep07(2017)150 (2017).
  * Maldacena and Qi [2018] J. Maldacena and X.-L. Qi, Eternal traversable wormhole, eprint (2018), arXiv:1804.00491 [hep-th] .
  * Cáceres  _et al._ [2021] E. Cáceres, A. Misobuchi, and R. Pimentel, Sparse SYK and traversable wormholes, Journal of High Energy Physics 2021, 11 (2021), arXiv:2108.08808 [hep-th] .
  * Nosaka and Numasawa [2021] T. Nosaka and T. Numasawa, Chaos exponents of SYK traversable wormholes, Journal of High Energy Physics 2021, 10.1007/jhep02(2021)150 (2021).
  * Nosaka and Numasawa [2023] T. Nosaka and T. Numasawa, On SYK traversable wormhole with imperfectly correlated disorders, Journal of High Energy Physics 2023, 10.1007/jhep04(2023)145 (2023).
  * Caldeira and Leggett [1981] A. O. Caldeira and A. J. Leggett, Influence of Dissipation on Quantum Tunneling in Macroscopic Systems, Phys. Rev. Lett. 46, 211 (1981).
  * Kamenev [2011] A. Kamenev,  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2011).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Reports on Progress in Physics 79, 096001 (2016).
  * Garc\́mathrm{i}a-Garc\́mathrm{i}a  _et al._ [2022c] A. M. Garc\́mathrm{i}a-Garc\́mathrm{i}a, Y. Jia, D. Rosa, and J. J. M. Verbaarschot, Replica Symmetry Breaking in Random Non-Hermitian Systems (2022c).
  * Kawabata  _et al._ [2023] K. Kawabata, A. Kulkarni, J. Li, T. Numasawa, and S. Ryu, Dynamical quantum phase transitions in Sachdev-Ye-Kitaev Lindbladians, Phys. Rev. B 108, 075110 (2023), arXiv:2210.04093 [cond-mat.stat-mech] .
  * Khramtsov and Lanina [2021] M. Khramtsov and E. Lanina, Spectral form factor in the double-scaled SYK model, Journal of High Energy Physics 2021, 10.1007/jhep03(2021)031 (2021).
  * Tarnopolsky [2019] G. Tarnopolsky, Large $q$ expansion in the Sachdev-Ye-Kitaev model, Phys. Rev. D 99, 026010 (2019), arXiv:1801.06871 [hep-th] .
  * Lin and Motrunich [2018] C.-J. Lin and O. I. Motrunich, Out-of-time-ordered correlators in a quantum Ising chain, Phys. Rev. B 97, 144304 (2018).
  * Craps  _et al._ [2020] B. Craps, M. De Clerck, D. Janssens, V. Luyten, and C. Rabideau, Lyapunov growth in quantum spin chains, Phys. Rev. B 101, 174313 (2020).
  * Jalabert  _et al._ [2018] R. A. Jalabert, I. Garc\́mathrm{i}a-Mata, and D. A. Wisniacki, Semiclassical theory of out-of-time-order correlators for low-dimensional classically chaotic systems, Phys. Rev. E 98, 062218 (2018).
  * Garcia-Mata  _et al._ [2023] I. Garcia-Mata, R. Jalabert, and D. Wisniacki, Out-of-time-order correlations and quantum chaos, Scholarpedia 18, 55237 (2023).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103 (2015), arXiv:1411.0660 [cond-mat.dis-nn] .
  * Maldacena  _et al._ [2020] J. Maldacena, G. J. Turiaci, and Z. Yang, Two dimensional nearly de Sitter gravity (2020), arXiv:1904.01911 [hep-th] .
  * Cotler  _et al._ [2020] J. Cotler, K. Jensen, and A. Maloney, Low-dimensional de Sitter quantum gravity, Journal of High Energy Physics 2020 (2020).
  * Aalsma and Shiu [2020] L. Aalsma and G. Shiu, Chaos and complementarity in de Sitter space, Journal of High Energy Physics 2020, 10.1007/jhep05(2020)152 (2020).
  * Gu and Kitaev [2019] Y. Gu and A. Kitaev, On the relation between the magnitude and exponent of OTOCs, Journal of High Energy Physics 2019, 10.1007/jhep02(2019)075 (2019).