# Operator growth in open quantum systems: lessons from the dissipative SYK

Budhaditya Bhattacharjee, Xiangyu Cao, Pratik Nandy, Tanay Pathak

► To cite this version:

Budhaditya Bhattacharjee, Xiangyu Cao, Pratik Nandy, Tanay Pathak. Operator growth in open quantum systems: lessons from the dissipative SYK. Journal of High Energy Physics, 2023, 03 (3), pp.054. 10.1007/JHEP03(2023)054. hal-03921151v2

HAL Id: hal-03921151

https://hal.science/hal-03921151v2

Submitted on 29 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Operator growth in open quantum systems: lessons from the dissipative SYK

Budhaditya Bhattacharjee,ª Xiangyu Cao,ᵇ,<fnref n="1" /> Pratik Nandyc,<fnref n="1" /> and Tanay Pathakª

**ABSTRACT:** We study the operator growth in open quantum systems with dephasing dissipation terms, extending the Krylov complexity formalism of [1]. Our results are based on the study of the dissipative $q$-body Sachdev-Ye-Kitaev (SYK$_q$) model, governed by the Markovian dynamics. We introduce a notion of “operator size concentration” which allows a diagrammatic and combinatorial proof of the asymptotic linear behavior of the two sets of Lanczos coefficients ($a_n$ and $b_n$) in the large $q$ limit. Our results corroborate with the semi-analytics in finite $q$ in the large $N$ limit, and the numerical Arnoldi iteration in finite $q$ and finite $N$ limit. As a result, Krylov complexity exhibits exponential growth following a saturation at a time that grows logarithmically with the inverse dissipation strength. The growth of complexity is suppressed compared to the closed system results, yet it upper bounds the growth of the normalized out-of-time-ordered correlator (OTOC). We provide a plausible explanation of the results from the dual gravitational side.

**KEYWORDS:** Field Theories in Lower Dimensions, Quantum Dissipative Systems, Holography and Condensed Matter Physics (AdS/CMT)

ARXIV ePRINT: 2212.06180

<footnote>ªCentre for High Energy Physics, Indian Institute of Science,
C.V. Raman Avenue, Bangalore 560012, India</footnote>

<footnote>ᵇLaboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS,
Sorbonne Université, Université de Paris,
Paris 75005, France</footnote>

<footnote>ᶜCenter for Gravitational Physics and Quantum Information,
Yukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan

E-mail: budhadityab@iisc.ac.in, xiangyu.cao@ens.fr,
pratik@yukawa.kyoto-u.ac.jp, tanaypathak@iisc.ac.in</footnote>

<footnote n="1">¹Corresponding author.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Krylov complexity in open systems: general approaches</td><td>3</td></tr><tr><td>2.1</td><td>Using the closed Krylov basis</td><td>4</td></tr><tr><td>2.2</td><td>Arnoldi iteration</td><td>4</td></tr><tr><td>2.3</td><td>The moment method</td><td>5</td></tr><tr><td>3</td><td>The dissipative SYK</td><td>7</td></tr><tr><td>4</td><td>Large q exact result</td><td>8</td></tr><tr><td>4.1</td><td>Operator size concentration</td><td>8</td></tr><tr><td>4.2</td><td>Consequence on moments</td><td>12</td></tr><tr><td>5</td><td>Finite q numerics</td><td>13</td></tr><tr><td>5.1</td><td>Large N: operator size distribution</td><td>13</td></tr><tr><td>5.2</td><td>Large N: Arnoldi iteration</td><td>15</td></tr><tr><td>5.3</td><td>Arnoldi iteration at finite N</td><td>16</td></tr><tr><td>6</td><td>Operator growth and Krylov complexity</td><td>17</td></tr><tr><td>6.1</td><td>General argument</td><td>18</td></tr><tr><td>6.2</td><td>Exact solution</td><td>19</td></tr><tr><td>6.3</td><td>Bound on chaos</td><td>22</td></tr><tr><td>6.4</td><td>Holographic interpretation</td><td>22</td></tr><tr><td>7</td><td>Conclusion</td><td>23</td></tr><tr><td></td><td>A Derivation of eq. (3.4)</td><td>24</td></tr></table>

## 1 Introduction

A proper understanding of operator growth in dissipative systems is of fundamental interest. This is because dissipative phenomena are ubiquitous in nature and for any practical purposes, it is desirable to understand the effect of the environment on the system itself. One can consider the system plus the environment as a whole system that follows the unitary dynamics [2]. Tracing out the environment leads the system dynamics to non-unitary. However, for the Markovian environment, the system dynamics can be efficiently formulated in terms of Lindbladian evolution [3, 4]. The Lindbladian consists of two parts, a part usually known as the Hermitian Liouvillian which governs the unitary closed system
dynamics, and another part which includes the information of the environment. This part
is not Hermitian and thus makes the whole evolution non-unitary.

There has been a growing interest to study operator growth in systems whose dynamics
are governed by such Lindbladian. Especially it is of utmost interest how the environment
affects the chaotic and integrable nature of the system [5–7]. Several recent studies have
taken promising approaches using different probes, namely, the spectral statistics [8–11],
operator-size distribution [12–14], the out-of-time-ordered correlator (OTOC) [15–18], and
Krylov (K-)complexity [1, 19, 20].

In this paper, we aim to study the operator growth in the dissipative Sachdev-Ye-Kitaev (SYK) model from the point of view of the K-complexity [1] (for an incomplete list of studies, see [21–38] and the references therein). The SYK model is a (0+1)-dimensional quantum mechanical model which consists of *N* fermions where *q* of them are interacting at a time [39]. The model is integrable for *q* = 2 and chaotic for *q* > 2. Moreover, the model is exactly solvable, and it is particularly analytically tractable in the large *N* and the large *q* limit [40]. In particular, at low temperatures, it shows an emergence of conformal symmetry which allows writing a Schwarzian action [40]. This action shares close similarity with the effective action in AdS₂ gravity [41], for which it received considerable attention from the holographic side [42–44]. The chaotic nature of this model is well captured by OTOC, satisfying the chaos bound [45]. Several recent studies have also studied K-complexity in SYK analytically for closed system [1, 33] and numerically for closed [21, 22] and open/dissipative systems [20]. The numerical results are performed for finite *N* and finite *q* (especially *q* = 4) [20–22], where the Krylov bases have been constructed using the Lanczos algorithm [46]. For the dissipative systems, one can implement the algorithms suitable to non-unitary dynamics, namely the Arnoldi iteration [47], which was studied in [19]. The algorithms are numerical, thus one needs to restrict to a finite number of fermions. However, to see the universal properties of the Lanczos coefficients, especially its asymptotic growth in the presence of the environment, analytical studies are highly desirable. Here we report some analytical results in the large *N* and the large *q* limit for the dissipative SYK which has not been reported before.

The dissipative model we choose is proposed by Kulkarni, Numasawa and Ryu [48] (see [49] for a similar model). Recently, it has been shown that the model is dual to a non-Hermitian two-site SYK model at low temperature, connected by a Keldysh wormhole [50]. It is closely analogous to a Maldacena-Qi [51] coupled SYK model, albeit with some significant differences. To establish our results, we take several different approaches. First, we give a diagrammatic and combinatorial proof of our claim eq. (4.2) and eq. (4.3) for the SYK model. The proof relies on the 1/q expansion, introducing the notion of “operator size concentration”. Then we give an analytical check by generalizing the method of moments [52]. Previously, a specific version of this method was applied to closed systems [1, 33]. We extract the analytical form of two sets of Lanczos coefficients, which correspond to the primary diagonal and primary off-diagonal elements of the Lindbladian matrix written in Krylov basis [19]. We find that only the diagonal elements of the Lindbladian are sensitive to dissipation. We confirm our analytical results by semi-analytics in finite q and large N limit, and by directly implementing the numerical Arnoldi iteration
in finite $q$ and finite $N$ limit. We emphasize that these three different approaches have their own advantages and disadvantages, in lieu of capturing the K-complexity behavior, for different types of systems. It is demonstrated that the large $q$ limit of the dissipative SYK model is a unique scenario where these approaches are equally viable.

The consistencies of all the results allow us to conjecture an operator growth hypothesis for the open systems, given by eq. (6.1). The analytical expressions of $a_n$ and $b_n$ are then used to derive analytical expressions for the Krylov basis wavefunctions, Krylov complexity, and Krylov variance. These expressions correspond to a special sub-class of the dissipative SYK model, where the dissipation strength is quadratically related to the random coupling $J$. For weak dissipation, this model is a good approximation to the one considered in the rest of the manuscript. As a result, the Krylov complexity exhibits suppressed growth and saturation compared to the exponential growth for the closed systems. We identify the saturation timescale which grows logarithmically to the inverse of the dissipation strength. Although our results are based on the large $q$ analysis of the SYK model, we believe the claim of eq. (6.1) to hold for any generic open quantum systems. In fact, some analysis on spin systems [19, 20] supports our claim. Based on the above observation, we argue the growth of K-complexity must upper bounds the growth of normalized out-of-time-ordered (OTOC), which might generalize the chaos bound [45] in open systems. Finally, we speculate a holographic interpretation of our results from the dual gravity side.

The manuscript is structured as follows. In section 2 we review a few recently explored approaches to the K-complexity in open systems. In section 3, we review the dissipative SYK model which will be our main playground. It was introduced in [48], and more details can be found on [50, 53]. In section 4, we derive the main analytical results in the large $q$ limit. Section 5 presents numerical evidence that these results still hold as an excellent approximation away from large $q$. In 6 we derive the consequences in terms of the main results in terms of K-complexity growth and discuss implications on scrambling. We also give a heuristic explanation of our results from the gravitational side. Finally, section 7 concludes the discussion with some open questions.

## 2 Krylov complexity in open systems: general approaches

The Krylov (K-)complexity of operators is mostly studied in closed Hamiltonian systems. Extending it to open systems is not straightforward at all, and so far there is no consensus on how one should proceed in general. Here, we restrict our attention to a restricted class of open dynamics that are described by a Lindblad master equation which admits the infinite temperature (maximally mixed) state $\rho_\infty \propto I$ as a stationary state and focus on the dynamics of operators in this ensemble. Even then, a few approaches are possible (and some have been explored), as we review in this section.

To set the common stage, let us recall the Lindblad equation in its general form [3, 4]:

$$ \dot{\rho} = -i[H, \rho] + \sum_k \left[ L_k \rho L_k^\dagger - \frac{1}{2} \{ L_k^\dagger L_k, \rho \} \right]. \qquad (2.1) $$
Here $H$ and $\rho$ are the Hamiltonian and the density matrix of the system, and $L_k$'s as the jump operators. As mentioned above, we shall assume that

$$ \sum_k \left[ L_k \rho_\infty L_k^\dagger - \frac{1}{2} \{ L_k^\dagger L_k, \rho_\infty \} \right] = 0, \qquad (2.2) $$

where $\rho_\infty \propto I$ represents the infinite-temperature state. In this stationary state, it is reasonable to consider the operator dynamics, governed by

$$ \mathcal{O}(t) = e^{i\mathcal{L}t}\mathcal{O}, \qquad (2.3) $$

where the Heisenberg-picture Lindbladian $\mathcal{L}$ is given by [19]

$$ \mathcal{L} = \mathcal{L}_H + \mathcal{L}_D, \quad \mathcal{L}_H\mathcal{O} = [H, \mathcal{O}], \quad \mathcal{L}_D\mathcal{O} = -i \sum_k \left[ \mp L_k^\dagger \mathcal{O} L_k - \frac{1}{2} \{ L_k^\dagger L_k, \mathcal{O} \} \right]. \quad (2.4) $$

Here, we view $\mathcal{L}$ as a superoperator acting on the space of operators (the minus sign has to be taken when both $L_k$ and $\mathcal{O}$ are fermionic operators [20]). It is written as a sum of a unitary Liouvillian $\mathcal{L}_H$ and a dissipation contribution $\mathcal{L}_D$. The space of operators is also endowed with the inner product

$$ (A|B) := \mathrm{Tr}[\rho_\infty A^\dagger B] = \mathrm{Tr}[A^\dagger B]/\mathrm{Tr}[1]. \qquad (2.5) $$

Then, the autocorrelation function of an operator $\mathcal{O}$ can be defined similarly as in the closed-system context:

$$ \mathcal{C}(t) = (\mathcal{O}|\mathcal{O}(t)) = (\mathcal{O}|e^{i\mathcal{L}t}|\mathcal{O}). \qquad (2.6) $$

Without loss of generality, we shall assume that $\mathcal{O}$ has unit norm: $(A|B) = 1$.

## 2.1 Using the closed Krylov basis

Ref. [8] proposed a simple way to extend the notion of K-complexity to the above open setting: we take the Krylov basis generated by the repeated action of the closed-system Liouvillian $\mathcal{L}_H$. Recall that this is an orthonormal basis $\{\mathcal{O}_0, \dots, \mathcal{O}_n, \dots\}$ that spans the Krylov subspaces:

$$ \mathrm{span}(\mathcal{O}_0, \dots, \mathcal{O}_n) = \mathrm{span}(\mathcal{O}, \mathcal{L}_H \mathcal{O}, \dots, \mathcal{L}_H^n \mathcal{O}). $$

Assuming that the open-system operator dynamics remain in such a Krylov subspace, one can define its K-complexity in exactly the same way. Now, on this basis the full Lindbladian is represented as a general matrix: while $\mathcal{L}_H$ becomes tridiagonal, $\mathcal{L}_D$ enjoys no special property in general. Nevertheless, the authors of [8] put forward conjectures on the matrix elements $(\mathcal{O}_n|\mathcal{L}_D|\mathcal{O}_m)$ motivated by numerical study, which we shall discuss below.

## 2.2 Arnoldi iteration

An alternative method, explored in ref. [19], considers another orthonormal basis $\{\mathcal{V}_0, \dots, \mathcal{V}_n, \dots\}$, which corresponds to the Krylov subspaces generated by the full open-system Lindbladian:

$$ \mathrm{span}(\mathcal{V}_0, \dots, \mathcal{V}_n) = \mathrm{span}(\mathcal{O}, \mathcal{L}\mathcal{O}, \dots, \mathcal{L}^n\mathcal{O}). $$
This basis, as well as the representation of the Lindbladian in it, can be generated by the Arnoldi iteration that we recall now.

We start with an initial normalized vector $V_0 \propto O$. For $k = 1, 2, \dots$, we construct

$$ |\mathcal{U}_k\rangle = \mathcal{L}|V_{k-1}\rangle. \qquad (2.7) $$

Then, for $j = 0$ to $n-1$, we perform the following iterations:

$$ 1. h_{j,k-1} = (V_j|\mathcal{U}_k). $$

$$ 2. |\tilde{\mathcal{U}}_k\rangle = |\mathcal{U}_k\rangle - \sum_{j=0}^{k-1} h_{j,k-1}|V_j\rangle. $$

$$ 3. h_{k,k-1} = \sqrt{(\tilde{\mathcal{U}}_k|\tilde{\mathcal{U}}_k)}. $$

(Stop if $h_{k,k-1} = 0$.) Otherwise define $V_k$ as

$$ |V_k\rangle = \frac{|\tilde{\mathcal{U}}_k\rangle}{h_{k,k-1}}. \qquad (2.8) $$

As a result, the Lindbladian is transformed into an upper Hessenberg form in the Krylov basis (sometimes called Arnoldi basis):

$$ \mathcal{L} \sim \begin{pmatrix} h_{0,0} & h_{0,1} & h_{0,2} & \cdots & \cdots & h_{0,n} \\ h_{1,0} & h_{1,1} & h_{1,2} & \cdots & \cdots & h_{1,n} \\ 0 & h_{2,1} & h_{2,2} & h_{2,3} & \cdots & \cdots \\ \cdots & 0 & h_{3,2} & \cdots & \cdots & h_{n-1,n} \\ 0 & \cdots & 0 & \cdots & \cdots & h_{n-1,n} \\ 0 & 0 & \cdots & 0 & h_{n,n-1} & h_{n,n} \end{pmatrix}, \qquad (2.9) $$

where $h_{m,n} = (V_m|\mathcal{L}|V_n)$. When $\mathcal{L}$ is a Hermitian superoperator, the Arnoldi iteration reduces to the Lanczos algorithm, and the above matrix becomes real symmetric and tridiagonal. However, for a general Lindbladian, these properties are not guaranteed. Thus, while the basis $V_k$ allows to define a version of K-complexity, the Lindbladian is still relatively involved.

## 2.3 The moment method

A distinct approach is to generalize the “moment method” [52], which was previously studied in closed SYK [1], and its higher-order corrections [33]. Rather than constructing a basis, this method focuses on the auto-correlation function, or its derivatives around $t = 0$, known as the moments:

$$ \mathcal{C}(t) \stackrel{t \ge 0}{=} \sum_{n=0}^{\infty} m_n \frac{(it)^n}{n!}. \qquad (2.10) $$

Now, assuming that the superoperator generating the dynamics is Hermitian, the Lanczos coefficients and the moments are related by a nonlinear transform. In terms of the Lanczos
coefficients, the moments, encoded in the Green function, can be obtained as a continued fraction expansion [1, 52]:

$$ \sum_{n=0}^{\infty} z^n m_n = \frac{1}{1 - a_0 z - \frac{b_1^2 z^2}{1 - a_1 z - \frac{b_2^2 z^2}{1 - a_2 z - \dots}}} \quad (2.11) $$

From the moments, the Lanczos coefficients can be calculated by a recursive algorithm [52]:

$$ \begin{align} M_k^{(0)} &= (-1)^k m_k, \quad L_k^{(0)} = (-1)^{k+1} m_{k+1}, \notag \\ M_k^{(n)} &= L_k^{(n-1)} - L_{n-1}^{(n-1)} \frac{M_{n-1}^{(n-1)}}{M_k^{(n-1)}}, \notag \\ L_k^{(n)} &= \frac{M_{k+1}^{(n)}}{M_n^{(n)}} - \frac{M_{n-1}^{(n-1)}}{M_k^{(n-1)}}, \quad k \ge n, \notag \\ b_n &= \sqrt{M_n^{(n)}}, \quad a_n = -L_n^{(n)}. \tag{2.12} \end{align} $$

Note that the above relations are a bit more general than those used in [1], which are valid for the case where the odd moments vanish (this is the case for closed systems and a Hermitian operator $O$, for example). Formally, we can apply the algorithm (2.12) to the moments of an open-system autocorrelation function, ignoring the above-mentioned assumption. Generically — that is, barring the occurrence of division-by-zero errors — we are guaranteed to obtain two sequences of complex numbers, $\{b_n\}_{n=1}^{\infty}$, $\{a_n\}_{n=0}^{\infty}$. However, it should be emphasized that they are *not* the output of a Lanczos algorithm. In fact, the moment method generates a *non-orthonormal* basis of the Krylov subspaces, in which $\mathcal{L}$ is represented as a symmetric tridiagonal matrix

$$ \mathcal{L} \sim \begin{pmatrix} a_0 & b_1 \\ b_1 & a_1 & b_2 \\ b_2 & a_2 & b_3 \\ & b_3 & \dots \\ & & \dots \\ & & & \dots \end{pmatrix}. \quad (2.13) $$

Thus, although we simplified the matrix representation of $\mathcal{L}$, it becomes problematic to define the notion of K-complexity using the non-orthonormal basis, which can be in practice quite singular (with respect to the inner product).

To conclude, we reviewed three different attempts to extend the K-complexity to an open-system context. None of them is entirely satisfactory; each of them solves certain problems by making others worse. There seems to be a list of “wishes” that cannot be simultaneously fulfilled in general. The main point of this paper is to show that all wishes can come true in *certain* systems.
### 3 The dissipative SYK

In this section, we review and introduce a concrete system, which is a dissipative version of the Sachdev-Ye-Kitaev (SYK) model [39, 42, 43]. The SYK model is a (0+1)-dimensional quantum mechanical model which consists of N fermions where q of them are interacting at a time. The SYK$_q$ Hamiltonian is given by [39]

$$H = i^{q/2} \sum_{i_1 < \cdots < i_q} J_{i_1 \cdots i_q} \psi_{i_1} \cdots \psi_{i_q}, \qquad (3.1)$$

where $\psi_i$'s are Majorana fermions satisfying $\{\psi_a, \psi_b\} = \delta_{ab}$, and $J_{i_1 \cdots i_q}$ are random couplings, drawn from some Gaussian ensemble with zero mean $\langle J_{i_1 \cdots i_q} \rangle = 0$ and variance $\langle J_{i_1}^2 \cdots J_{i_q}^2 \rangle = (q-1)!J^2/N^{q-1}$. The constant $J$ is a dimensionful parameter and sets the energy scale of the Hamiltonian. We will be particularly interested in the limit where we send $N \to \infty$, and then $q \to \infty$. In that limit, the model is quantum chaotic yet enjoys a high degree of analytical tractability; numerical studies are required more or less away from this limit. When $q$ is large, it is convenient to define a rescaled coupling constant:

$$\mathcal{J}^2 = 2^{1-q} q J^2. \qquad (3.2)$$

To introduce openness, we consider the following jump operators [48]

$$L_i = \sqrt{\mu} \psi_i, \quad i = 1, \cdots, N. \qquad (3.3)$$

where $\mu \ge 0$ is the dissipation strength or that of the system-environment coupling. The dissipation operators are fermionic and Hermitian. The Hermiticity ensures that the infinite-temperature density matrix $\rho_\infty \propto I$ [48–50] is a stationary state. It will be useful to note that, the action of the dissipative part of the Lindbladian $\mathcal{L}_D$ is very simple on “Majorana string” operators. For any $s > 0$ and $1 \le i_1 < \cdots < i_s \le N$, we have

$$\mathcal{L}_D(\psi_{i_1} \cdots \psi_{i_s}) = i\mu s (\psi_{i_1} \cdots \psi_{i_s}). \qquad (3.4)$$

See appendix A for the derivation. In other words, the effect of the dissipation is to annihilate a Majorana string with a rate proportional to its size $s$ and to $\mu$ (recall that the operator evolution is generated by $i\mathcal{L}$). Therefore, $\mathcal{L}_D$ and $\mathcal{L}_H$ have competing effects as the latter makes the operator size grow. In the large $q$ limit, we rescale $\mu$ as well, so that

$$\tilde{\mu} = \mu q, \qquad (3.5)$$

is of order unity [48].

Often, the dynamics of the operators (as well as of the density matrix) is seen as taking place in the double Hilbert space $\mathcal{H} \otimes \mathcal{H}^* = \mathcal{H}_L \otimes \mathcal{H}_R$ through the lens of Choi-Jamiolkowski (CJ) isomorphism [54, 55]. In this regard, the operator dynamics is generated by the non-Hermitian Hamiltonian [50]:

$$H = iH_L^{\text{SYK}} - i(-1)^{q/2} H_R^{\text{SYK}} + i\mu \sum_j \psi_L^j \psi_R^j. \qquad (3.6)$$
This can be understood as two SYK (with relative phase between them) interacting with the coupling $H_{\text{int}} = i\mu \sum_j \psi_L^j \psi_R^j$. The state $|\mathcal{I}\rangle$, which is the image of $\rho_\infty = I$ under the CJ-isomorphism, is the ground state of $H_{\text{int}}$ and well as $H$. Note that this is slightly different than ref. [51], where the Hamiltonian is Hermitian (the imaginary factor is absent in the first two terms). Hence, the dissipative dynamics are equivalent to non-Hermitian dynamics governed by the above Hamiltonian.

## 4 Large $q$ exact result

We are in a position to state the main result of this work: for the dissipative SYK in the large $q$ limit with $O \propto \psi_1$, the three approaches in section 2 *exactly coincide*, in the sense that:

1. The closed-system Krylov basis $\{O_n\}$ is identical to the open-system one $\{V_n\}$.

2. The Arnoldi iteration yields a symmetric tridiagonal matrix, whose matrix elements are equal to the Lanczos coefficients obtained by the moment method. In particular, we have

$$h_{n,n-1} = h_{n-1,n} = b_n, \quad h_{n,n} = a_n. \qquad (4.1)$$

Furthermore, these coefficients are exactly calculated:

$$a_n = i\tilde{\mu}n + O(1/q), \quad \tilde{\mu} := \mu q, \qquad (4.2)$$

$$b_n = \begin{cases} \mathcal{J}\sqrt{2/q} & n=1, \\ \mathcal{J}\sqrt{n(n-1)} + O(1/q) & n>1. \end{cases} \qquad (4.3)$$

Here, the results on $b_n$ are identical to the closed-system one obtained in [1], which is expected from 1 above. The result on $a_n$ is new and shows that the diagonal elements are imaginary, and grow linearly with $n$ with a slope set by the dissipation strength [20].

In the rest of this section, we first present a direct derivation of the above results. Then, a nontrivial check of the Lanczos coefficients formulae (4.2) and (4.3) via the moment method will be exposed. Finally, numerical results away from the ideal limit will be discussed.

### 4.1 Operator size concentration

It turns out that the above statements are all corollaries of the following “operator size concentration” property of the Krylov basis elements of the closed SYK model in the large $q$ limit: the $n$-th basis operator $O_n$ is a linear combination of Majorana strings of the same size

$$O_n = \sum_{i_1<\dots<i_s} c_{i_1,\dots,i_s} \psi_{i_1} \dots \psi_{i_s} + O(1/q), \qquad (4.4)$$

where

$$s = n(q-2) + 1. \qquad (4.5)$$
Here $s$ is the size and $n$ is the generation, according to the nomenclature of ref. [56]. Indeed, recall (3.4) that Majorana strings of size $s$ are degenerate eigen-operators of the dissipative part of Lindbladian, with eigenvalue $i\mu s$. Thus, (4.4) implies

$$ \mathcal{L}_D \mathcal{O}_n = i\mu s \mathcal{O}_n = i\tilde{\mu}n \mathcal{O}_n + O(1/q), \qquad (4.6) $$

where $\tilde{\mu}$ is defined through (3.5). This, combined with the known action of $\mathcal{L}_H$, implies the statements above as well as (4.2) and (4.3) immediately.

It remains to prove operator size concentration (4.4), which concerns the closed SYK model in the large $q$ limit. It is known [1, 56] that one can study the operator dynamics and implement the Lanczos algorithm directly in the large $N$ limit using “open” melon diagrams. Our proof will rely on this combinatorial approach, which we briefly review. For simplicity we set $\mathcal{J} = 1/\sqrt{2}$ below.

Let us iteratively apply $\mathcal{L}_H$ to $\psi_1$ and describe graphically the operators generated. It will be convenient to decompose the Liouvillian as [27]

$$ \mathcal{L}_H = \mathcal{L}_+ + \mathcal{L}_-, \qquad (4.7) $$

where $\mathcal{L}_\pm$ are the operator-size increasing (decreasing, respectively) contribution to $\mathcal{L}_H$. In other words, $\mathcal{L}_+$ is defined such that for two Majorana strings of sizes $s$ and $t$,

$$ (\psi_{i_1} \dots \psi_{i_t} | \mathcal{L}_+ | \psi_{i_1} \dots \psi_{i_s}) = \begin{cases} (\psi_{i_1} \dots \psi_{i_t} | \mathcal{L}_H | \psi_{i_1} \dots \psi_{i_s}) & t > s \\ 0 & \text{otherwise.} \end{cases} \qquad (4.8) $$

$\mathcal{L}_-$ is then defined as the Hermitian conjugate of $\mathcal{L}_+$. By applying $\mathcal{L}_+$ instead of $\mathcal{L}_H$, we can focus on the new operators generated at each step.

At the first step,

$$ \mathcal{L}_+ \psi_1 \propto \begin{array}{c} \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) \\ \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) \end{array} = \begin{array}{c} \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) \\ \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) \end{array}. \qquad (4.9) $$

Here, we first represent $\mathcal{L}_H \psi_1$, which is a $(q-1)$-body operator generated from the $\psi_1$, by the half of a melon diagram as in [56]. Since $q$ is large, we omit all (grey) propagators but two representatives, and also the tip of the melon. Hence the diagram simplifies to a simple arc.

With this in mind, the next few operators generated will be the following:

$$ \mathcal{L}_+^2 \psi_1 \propto \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right), \qquad (4.10) $$

$$ \mathcal{L}_+^3 \psi_1 = c_3 \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) + c_4 \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right), \qquad (4.11) $$

$$ \mathcal{L}_+^4 \psi_1 = c_5 \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) + c_6 \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) + c_7 \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right) + c_8 \left( \begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_n \end{array} \right), \qquad (4.12) $$
and so on. Every time, the $\mathcal{L}_+$ acts on one of the Majorana's inside an arc, and transforms it into $(q-1)$ Majorana's represented by a new arc, which we view as a "child" of the former. We can equivalently view the operators appearing in $\mathcal{L}_+^n \psi_1$ as rooted and unmarked trees with $n$ vertices (a vertex = an arc). We remark that equations (4.10) and (4.12) are operator identities that are valid before averaging over the disordered couplings $J_{i_1 \dots i_q}$, modulo terms that will have negligible contributions to any disordered-averaged observables in the large $N$ limit. Each arc involves a dangling disorder line. The subleading terms are neglected by considering only open melon diagrams. The observables can be constructed by closing the open melon diagrams, that is, taking the inner product with another operator and averaging over disorder.

So far, we have been elusive about the prefactors. A neat way to keep track of them is to order the arcs in the diagrams, such that a vertex comes always before its children. For example, the diagram with prefactor $c_7$ above can be ordered in three ways:

$$c_7 = \begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
    \node (n1) at (0,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}4 \\ 3 \end{array}$};
    \end{tikzpicture}$};
    \node (n2) at (1,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}4 \\ 2 \end{array}$};
    \end{tikzpicture}$};
    \node (n3) at (2,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}3 \\ 2 \end{array}$};
    \end{tikzpicture}$};
    \draw (n1) -- (n2);
    \draw (n2) -- (n3);
\end{tikzpicture}.$$

In other words, the amplitude of an unmarked diagram in $\mathcal{L}_+^n \psi_1$ is given by the ways it can be built by adding one arc at a time. Note that the terms of the right-hand side are the same operator; the different orderings enumerate its amplitude (multiplicity) in $\mathcal{L}_+^n \psi_1$. As a consequence,

$$\mathcal{L}_+^n \psi_1 = \sum \text{[ordered diagrams of } n \text{ vertices]} . \qquad (4.13)$$

We now come to the essential claim of this section. For any $n \ge 1$, we have

$$\mathcal{L}_{-} \mathcal{L}_{+}^{n+1} \psi_1 = \frac{1}{2} n(n+1) \mathcal{L}_{+}^{n} \psi_1 . \qquad (4.14)$$

This identity can be proved diagrammatically. By (4.13) and the definition of $\mathcal{L}_-$, the left-hand side is a sum of ordered diagrams of $n+1$ arcs with one childless arc marked. For example,

$$\mathcal{L}_{-} \mathcal{L}_{+}^{3+1} \psi_1 = \mathcal{L}_{-} \begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
    \node (n1) at (0,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}4 \\ 3 \end{array}$};
    \end{tikzpicture}$};
    \node (n2) at (1,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}4 \\ 2 \end{array}$};
    \end{tikzpicture}$};
    \node (n3) at (2,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}3 \\ 2 \end{array}$};
    \end{tikzpicture}$};
    \draw (n1) -- (n2);
    \draw (n2) -- (n3);
\end{tikzpicture} + \dots = \begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
    \node (n1) at (0,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}4 \\ 3 \end{array}$};
    \end{tikzpicture}$};
    \node (n2) at (1,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}4 \\ 2 \end{array}$};
    \end{tikzpicture}$};
    \node (n3) at (2,0) {$\begin{tikzpicture}[baseline={([yshift=-0.5ex]current bounding box.center)}]
        \draw (0,0) circle (0.5);
        \node at (0,0) {$\begin{array}{c}3 \\ 2 \end{array}$};
    \end{tikzpicture}$};
    \draw (n1) -- (n2);
    \draw (n2) -- (n3);
\end{tikzpicture} + \dots . \qquad (4.15)$$

In the second line, the marked childless arc is represented by a red dashed line. One should view it as being removed by the action of $\mathcal{L}_-$.

Now, to prove (4.14), we shall show that the removal of the marked arc gives rise to an $n(n+1)/2$-to-one correspondence between the ensemble of marked and ordered diagrams with $n+1$ arcs, and that of ordered diagrams with $n$ arcs. To do this, consider a marked and ordered diagram of $3+1$ arcs, say,

<figure><img src="image_2.png" /></figure>
following ordered diagram with 3 arcs:

<figure><img src="image_3.png" /></figure>

This removal map is not one-to-one, since we lost information about the marked arc, namely, its parent's index $p \in \{1, \dots, n\}$, and its order $q \in \{p+1, \dots, n+1\}$ (since it must come later than its parent). In the above example, $(p, q) = (1, 2)$ in the example. It is not hard to see that the datum of

$$ (p, q) \text{ s.t. } p \in \{1, \dots, n\}, q \in \{p+1, \dots, n+1\}, \qquad (4.16) $$

allows to reconstruct uniquely the marked ordered diagram with $n+1$ arcs from any ordered diagram of $n$ arcs. For example,

<figure><img src="image_4.png" /></figure>

Observe that removing the marked arc on the right-hand side gives back the diagram on the left-hand side. Since there are

$$ 1 + 2 + \dots + n = \frac{1}{2}n(n+1) $$

such tuples satisfying (4.16), the removal map is a $n(n+1)/2$-to-one correspondence. This concludes the combinatorial proof of (4.14).

The key consequence of (4.14) is that, the Krylov basis generated by $\mathcal{L}_H$ and $\psi_1$ is essentially $\{\mathcal{L}_+^n \psi_1\}$, up to normalization:

$$ \mathcal{O}_n \propto \mathcal{L}_+^n \psi_1. \qquad (4.17) $$

To see this, we first notice that $\{\mathcal{L}_+^n \psi_1\}$ is an orthogonal set, since $\mathcal{L}_+^n \psi_1$ is a linear combination of Majorana strings of length $s = n(q-2)+1$, and Majorana strings of distinct lengths are orthogonal under the inner product (2.5). Moreover, by (4.7) and (4.14), for $n > 1$,

$$ \mathcal{L}_H \mathcal{L}_+^n \psi_1 = (\mathcal{L}_+ + \mathcal{L}_-) \mathcal{L}_+^n \psi_1 = \mathcal{L}_+^{n+1} \psi_1 + \frac{1}{2} n(n-1) \mathcal{L}_+^{n-1} \psi_1. \qquad (4.18) $$

The edge cases $n=0,1$ can be taken care of by explicit calculation:

$$ \mathcal{L}_H \psi_1 = \mathcal{L}_+ \psi_1, \qquad \mathcal{L}_H \mathcal{L}_+ \psi_1 = \mathcal{L}_+^2 \psi_1 + \frac{1}{q} \psi_1. \qquad (4.19) $$

Therefore, the Liouvillian $\mathcal{L}_H$ is represented as a tridiagonal matrix in the basis $\{\mathcal{L}_+^n \psi_1\}$. Combined with its orthogonality, we conclude that it must be identical to the Krylov
basis up to normalization, (4.17). In fact, (4.18) implies that the Lanczos coefficients $b_n = \sqrt{n(n-1)/2}$ for $n > 1$, as announced in (4.3) (recall that we set $\mathcal{J} = 1/\sqrt{2}$).<fnref n="1" /> Thus, we have provided an independent, diagrammatic derivation of this nontrivial result (the derivation in [1] relied on the moment method and knowledge about tangent numbers, see also below).

Eq. (4.17) implies operator size concentration (4.4) immediately: since $\mathcal{L}_+^n \psi_1$ is by construction a linear combination of Majorana strings of length $s = n(q-2)+1$, the same is true for the $n$-th Krylov basis element of the large $q$ SYK model. This concludes the derivation of the main results announced at the beginning of the section.

## 4.2 Consequence on moments

A consequence of the main results is that the Lanczos coefficients (4.2) and (4.3) are related to the autocorrelation function of the dissipative SYK in the large $q$ limit. The latter has been computed explicitly by solving a Schwinger-Dyson equation [48]:<fnref n="2" />

$$ \mathcal{C}(t) = 1 + \frac{1}{q}g(t) + O(1/q^2), \qquad (4.20) $$

$$ g(t) = \log \left[ \frac{\alpha^2}{\mathcal{J}^2 \cosh^2(\alpha t + \gamma)} \right], \quad t > 0, \qquad (4.21) $$

where we recall that $\mathcal{J}^2 = 2^{1-q}q\mathcal{J}^2$, $\tilde{\mu} = \mu q$, and the other variables are given by

$$ \alpha = \sqrt{\left(\frac{\tilde{\mu}}{2}\right)^2 + \mathcal{J}^2}, \quad \gamma = \operatorname{arcsinh}\left(\frac{\tilde{\mu}}{2\mathcal{J}}\right). \qquad (4.22) $$

In the closed-system case $\mu = 0$, we recover the known result $g(t) = 2 \ln \operatorname{sech}(\mathcal{J}t)$ [40], and the associated moments are essentially the tangent numbers, and the Green function is known to have a continued fraction expansion given by the $b_n$'s [1]. When $\mu > 0$, the moments are considerably more involved. Setting $\mathcal{J} = 1$ without loss of generality, the moments have the large $q$ expansion

$$ m_n = \frac{2}{q} \tilde{m}_n + O(1/q^2), \quad n \ge 1, \qquad (4.23) $$

and $\tilde{m}_n$ is a polynomial of $u := i\tilde{\mu}$. For example, $\tilde{m}_1 = u/2$,

$$ \tilde{m}_2 = 1, \qquad (4.24) $$

$$ \tilde{m}_3 = u, \qquad (4.25) $$

$$ \tilde{m}_4 = u^2 + 2, \qquad (4.26) $$

<footnote n="1">¹The square-root comes from the appropriate normalization of $\mathcal{L}_+^n \psi_1$, to get $\mathcal{O}_n$. The way to proceed with this is to consider the norm, $\|\mathcal{L}_+^n \psi_1\|^2 = (\psi_1|\mathcal{L}_-^n \mathcal{L}_+^n|\psi_1) = \frac{n(n-1)}{2} \frac{(n-1)(n-2)}{2} \dots (\psi_1|\psi_1)$, which gives the normalization squared.</footnote>

<footnote n="2">²We remark that the large $q$ expansion above is valid for $t \ll t_c$, where $t_c$ is such that $|g(t_c)| = O(q)$, or $t_c \sim q$, for large $q$. Beyond that time scale, the expansion (4.20) breaks down, and one would need a different approximation [53].</footnote>
$$ \tilde{m}_5 = u^3 + 8u, \qquad (4.27) $$

$$ \tilde{m}_6 = u^4 + 22u^2 + 16, \qquad (4.28) $$

$$ \tilde{m}_7 = u^5 + 52u^2 + 136, \qquad (4.29) $$

$$ \tilde{m}_8 = u^6 + 114u^4 + 720u^2 + 272, \qquad (4.30) $$

and so on. The coefficients appear to coincide with the so-called triangle “$T(n, k)$”, see [57] and references therein for further mathematical facts about it. Our main result then implies the following continued fraction formula

$$ \sum_{n \ge 0} \tilde{m}_{n+2z^n} = \frac{1}{1 - uz - \frac{1 \times 2z^2}{1 - 2uz - \frac{2 \times 3z^2}{1 - 3uz - \frac{3 \times 4z^2}{1 - 4uz - \dots}}} \qquad (4.31) $$

One may check this formula order by order in $z$. We are not aware of an independent “elementary” proof. See however discussion around eq. (6.19) below.

# 5 Finite q numerics

In this section, we turn to examine how the simple “ideal” scenario established above at the $q \to \infty$ limit is affected at finite $q$.

## 5.1 Large N: operator size distribution

Our main results at the large $q$ limit rely on the operator size concentration property of the closed-system Krylov basis. A natural way to anticipate how these results survive at finite $q$ is to examine the operator size distribution [56, 58]. This can still be done in the large $N$ limit, by extending the diagrammatic approach to large $N$ operators [1, 56] described in the previous section to an arbitrary value of $q$: a finite $q$ imposes that every arc can have at most $q-1$ (direct) children. Coding this in a computer, we may generate all the large $N$ operators up to a certain size, which allows us to implement the Lanczos algorithm in the large $N$ limit and construct explicitly the Krylov basis operators, up to a certain value of $n$. (As the number of operators increases exponentially in $n$, we are limited to $n \approx 20$; fortunately, this is enough to observe the asymptotic behaviors.)

We then measure the operator size distribution of the Krylov basis operators $O_n$. In a finite size system, this is defined as follows [56]: writing $O_n$ as a sum of Majorana strings

$$ O_n = \sum_s \sum_{i_1 < \dots < i_s} c_{i_1, \dots, i_s} \psi_{i_1} \dots \psi_{i_s}, \qquad (5.1) $$

the size distribution is then given by

$$ P_s = \sum_{i_1 < \dots < i_s} |c_{i_1, \dots, i_s}|^2. \qquad (5.2) $$
<figure><img src="image_5.png" /><figcaption>**Figure 1.** (a) Operator size distribution of the Krylov basis operators $O_n$ in large $N$, $q = 4$ SYK model. In the main plot, we plot the standard deviation of the operator size divided by the average size of $O_n$ as a function of $n$. The inset shows the whole operator size distribution for a few Krylov basis operators (the value of $n$ is given by the color code which is the same as in the main). The dashed lines are a guide to the eye. The distribution is strongly peaked at the largest operator size $s = (q-2)n+1$ present in $O_n$, so that the standard deviation is $\le 0.016$ times the mean for all $n \le 17$, and appears to tend to a constant as $n$ increases. (b) Magnitude of the matrix elements of the matrix $h_{m,n}$ resulting from the Arnoldi iteration in the large $N$, $q = 4$ dissipative SYK model with $\tilde{\mu} = 1$ and $J = 1$. The presence of the secondary off-diagonal elements in the upper Hessenberg form indicates the finite $q$ effect.</figcaption></figure>

In the large $N$ representation, an operator is a sum of diagrams, each representing a sum of Majorana strings with definite size $s$, related to the number of arcs $n$ as $s = (q-2)n+1$. The weights $|c_{i_1, \dots, i_s}|^2$ depends on the random coupling coefficients $J_{i_1 \dots i_q}$. Yet, one can show that each diagram after disorder averaging and summing over fermion flavors, each diagram gives the same weight. Hence, the size distribution can be readily measured at large $N$. Using the size distribution, we can obtain readily the average operator size $\langle s \rangle := \sum_s s P_s$ as well as higher moments.

In figure 1(a), we plot the results for $q = 4$. We observe that the operator size concentration property holds exactly if and only if $n = 0, \dots, 4$: only the first five Krylov basis operators are a sum of Majorana's of a definite size $s = (q-2)n+1$. For larger $n$, the exact operator size concentration is indeed only a property of the large $q$ limit.<fnref n="3" /> Nevertheless, the distribution of operator size is highly peaked at $s = (q-2)n+1$, which is the maximal operator size present in $O_n$. The weight of operator sizes $s' < s$ decays exponentially fast in $s-s'$, see the inset of figure 1(a). As a result, the mean operator size $\langle s \rangle \approx (q-2)n+1$ grows linearly in $n$. Its standard deviation also grows linearly with $n$, yet much more slowly than the mean value. Our numerics indicates the standard deviation-mean ratio tends to be a small constant:

$$ \frac{\sqrt{\langle s^2 \rangle - \langle s \rangle^2}}{\langle s \rangle} \rightarrow C \approx 0.015, \qquad (5.3) $$

for $q = 4$ SYK. Data for other values of $q$ suggest that the constant decays as $C \approx 0.06/q$.

<footnote n="3">³We have checked a few other values of $q$ and observed that the property always breaks down at $n = 5$.</footnote>
In summary, we provided significant evidence indicating that operator size concentration is an excellent approximation in the SYK model for any values of $q$. We may further rationalize this optimism by observing that the operator size concentration holds exactly for $q = 2$ (since only Majorana strings of length 1 are generated) as well as for $q = \infty$; hence, it cannot be too wrong at intermediate values of $q$.

## 5.2 Large N: Arnoldi iteration

As a consequence of operator size concentration, we expect that the Arnoldi iteration should result in a matrix $h_{m,n}$ close to the “ideal scenario” [19, 20]. That is, to a good approximation,

1. $[h_{mn}]$ should be a symmetric tridiagonal matrix,

2. its diagonal elements are imaginary and grow linearly in $n$ with a pre-factor proportional to $\mu$:

$$h_{n,n} \approx i\mu\chi n. \qquad (5.4)$$

Here $\chi$ is some $\mu$-independent proportionality constant.

3. Its primary off-diagonal elements grow linearly in $n$:

$$h_{n-1,n} \approx \alpha n, \qquad (5.5)$$

where $\alpha$ is independent of $\mu$.

This is indeed what we found from our large $N$ numerics. In figure 1(b) we plot the magnitude of the matrix elements $h_{m,n}$ obtained from the Arnoldi iteration. As expected, the largest matrix elements are on the diagonal and primary off-diagonals; the other nonzero elements, that is, $h_{m,n}$ for $m < n - 1$, are about two orders of magnitudes smaller (both the dissipation rate and the unitary coupling constant is of order one, $\tilde{\mu} = J = 1$).

Next, let us examine more quantitatively how close the matrix is to the “ideal scenario”. In figure 2(a), we plot the primary off-diagonal elements $h_{n+1,n}$ (which are real), the imaginary part of the diagonal elements $\text{Im}(h_{n,n})$ (they turn out purely imaginary), and also an “error” defined as

$$\varepsilon_n^2 := |h_{n-1,n} - h_{n,n-1}|^2 + \sum_{k<n-1} |h_{k,n}|^2, \quad n > 0. \qquad (5.6)$$

By definition, the error vanishes if and only if $h$ is symmetric. We observe almost perfect linear growth in $n$ for both $h_{n+1,n}$ and $\text{Im}(h_{n,n})$. The growth rate of $\text{Im}(h_{n,n})$ is proportional to the dissipation rate $\tilde{\mu}$ while that of $h_{n+1,n}$ is almost independent of it (as long as $\tilde{\mu} \lesssim 1$). For the values of $n$ accessible to us, the error $\varepsilon_n$ is much smaller for order-unity dissipation $\tilde{\mu} = 1$. We also observe a rather fast growth of $\varepsilon_n$ as $n$ increases, and cannot rule out the possibility that the matrix $h_{m,n}$ deviates significantly from the ideal scenario for $n$ sufficiently large. However, for weaker dissipation, the errors appear to be well described by $\varepsilon_n \propto \mu n$, with a proportionality constant of order $10^{-2}$. Thus, the matrix should be close to the ideal scenario for $n \lesssim 10^2/\mu$. As we shall see below, the operator dynamics
<figure><img src="image_6.png" /><figcaption>Figure 2. (a) The diagonal and primary off-diagonal elements in the matrix $h_{m,n}$ resulting from the Arnoldi iteration in the large $N$, $q = 4$ dissipative SYK model with $\tilde{\mu} = 1, 0.1, 0.01$ and $J = 1$. The markers indicate the value of $\tilde{\mu}$. We also plot the quantity $\epsilon_n$ defined in (5.6) which measures the distance to the ideal scenario (it is multiplied by 100 for visibility). The diagonal elements and $\epsilon_n$ are rescaled by $\tilde{\mu}$ to display the collapse. (b) SYK$_4$ (i.e., $q = 4$) Lindbladian matrix (single realization) in Krylov (Arnoldi) basis for dissipation $\mu = 0.01$. We fix $J = 1$ and the system size $N = 18$. The presence of the secondary off-diagonal elements in the upper Hessenberg form is due to both $q$ and $N$ being finite, and to the absence of disorder averaging; compare with figure 1(b).</figcaption></figure>

will be confined in the region $n \lesssim 1/\mu$ under the ideal scenario, so the deviation from the ideal scenario at larger $n$ is irrelevant.

To conclude, the large $N$ numerics confirms that operator size concentration is a good approximation at finite $q$, and as a consequence, the Arnoldi iteration results in a matrix close to the “ideal scenario”: symmetric and tridiagonal, with diagonals $h_{n,n} \sim i\chi\mu n$ and primary off-diagonals $h_{n+1,n} \sim \alpha n$.

## 5.3 Arnoldi iteration at finite $N$

We now consider Arnoldi iteration for single realizations of the $q = 4$ SYK model at finite $N$, up to $N = 18$, starting from the (normalized) initial operator $\sqrt{2}\psi_1$. A sample of the resulting matrix $[h_{m,n}]$ is plotted in figure 2(b). It is still close to the ideal scenario, modulo the saturation effect induced by the finite system size (see below). The diagonal elements $h_{n,n}$ are imaginary and grow linearly with $n$ up to the finite size saturation. The growth rate compares well with the prediction (see also figure 3(a)):

$$h_{n,n} \approx i\mu s = i\mu(2n + 1). \qquad (5.7)$$

In other words, $\chi = q - 2 = 2$ in eq. (5.4). To obtain (5.7), we assume that operator size distribution is still a good approximation. Hence, we have $V_n \approx O_n$, and the operator size distribution of $V_n$ is concentrated at $(q-2)n+1 = 2n+1$. Then (5.7) follows from recalling how the dissipative Lindbladian acts on Majorana strings (3.4). The primary off-diagonals are also not exactly equal, $h_{n,n-1} \neq h_{n-1,n}$, see figure 3(b), yet the relative difference is small. The magnitudes of the other off-diagonal elements are also small compared to diagonal and primary off-diagonal elements which are manifested in figure 2(b).
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Behavior of the (a) primary diagonal and (b) off-diagonal elements respectively, for two different dissipation $\mu = 0.01$ and $\mu = 0.02$. The dotted lines in (a) are the fitted straight lines (5.7). The inset in (b) shows the difference $|h_{n,n-1} - h_{n-1,n}|$ grows with the dissipation. (c) The behavior of the diagonal elements (magnitude) for different system sizes. The black dotted line denotes the straight line (5.7) (with an offset introduced for visualization). (d) The saturation value (averaged for $n = 10$ to $n = 40$) of the diagonal elements for different system sizes. In all cases, our system size is $N = 18$ (single realization) and we choose $q = 4$, $J = 1$, and the fixed dissipation $\mu = 0.01$ in (c) and (d).</figcaption></figure>

Finally, in figure 3(c)–3(d), we exhibit the system size ($N$) dependence of the diagonal elements. In figure 3(c), we see a saturation plateau after an initial linear growth for $n \lesssim N/q$. The saturation value increases linearly with $N$, see figure 3(d). The finite size saturation happens since the Majorana string lengths are upper bounded by $N$, by (3.4). A similar behavior applies to the primary off-diagonal entries, as in the closed case.

## 6 Operator growth and Krylov complexity

In the preceding sections, we showed that — exactly at the large $N$ and large $q$ limit, to a good approximation beyond — the dissipative SYK model displays operator size concentration. As a result, the three approaches of section 2 coincide, and the Lindbladian acts on the Krylov basis, as a symmetric tridiagonal matrix whose diagonal and primary
off-diagonal elements have the following asymptotic growth:

$$a_n \approx i\chi\mu n, \quad b_n \approx \alpha n, \qquad (6.1)$$

This is also referred to as the “ideal scenario”. The coefficient $\alpha$ is the same as in the closed system, whereas $\chi$ depends on how the size of the Krylov basis operator $O_n$ grows with $n$. It may be a bit hasty to discuss how general it is. Even if we restrict to situations where the ideal scenario entails operator size concentration, it is in general unknown how the latter is respected in models away from large $N$; however, there is numerical evidence indicating that the ideal scenario remains a good approximation in quantum spin chains [19, 20]. We leave a detailed study to the future.

Nevertheless, we can readily examine the consequence of the ideal scenario on the K-complexity growth. For this, we expand the operator on the closed-system Krylov basis:

$$|O(t)) = \sum_n i^n \varphi_n(t) |O_n). \qquad (6.2)$$

Note that the norm of the operator is not conserved [20]: $\sum_n |\varphi_n(t)|^2 \neq 1$. Thus, we shall define the K-complexity as the average position $n$ of the normalized wavefunction:

$$K(t) = \frac{1}{\mathcal{Z}} \sum_n n |\varphi_n|^2, \quad \mathcal{Z} := \sum_n |\varphi_n|^2. \qquad (6.3)$$

Under the ideal scenario, the amplitudes satisfy the following differential equation

$$\partial_t \varphi_n(t) = ia_n\varphi_n(t) - b_{n+1}\varphi_{n+1}(t) + b_n\varphi_{n-1}(t). \qquad (6.4)$$

where $a_n$ and $b_n$ satisfy (6.1), and initially, the wavefunction is localized at the origin: $\varphi_n(t = 0) = \delta_{n,0}$.

## 6.1 General argument

To understand the asymptotic quantitative behavior in a simple way, we can make the Ansatz that $\varphi_n = \varphi(n)$ depends on $n$ smoothly, and approximate (6.4) by a PDE:

$$\partial_t \varphi + n(\chi\mu\varphi + 2\alpha\partial_n\varphi) = 0. \qquad (6.5)$$

This equation has a stationary ($t \to \infty$) solution:

$$\varphi_*(n) \propto e^{-n/\xi}, \quad \xi := \frac{2\alpha}{\chi\mu}. \qquad (6.6)$$

It has an exponential tail of width $\xi$, which is inverse proportional to the dissipation strength $\mu$, and proportional to the K-complexity growth rate in the closed system $2\alpha$. The continuum approach here is quantitatively correct if $\xi \gg 1$, i.e., in the weak dissipation regime. In this regime, the dissipation term is negligible for $n \ll \xi$. Now, initially, the wavefunction is localized at origin, so it will evolve as in the $\mu = 0$ case: it is known [1] that it will spread exponentially fast, and its average position $K(t) \sim e^{2\alpha t}$. This is true until the wavefunction reaches $K(t) \sim \xi$, where the dissipation term can no longer be neglected. In
fact, it suppresses any further growth of the K-complexity, and the wavefunction converges to the stationary solution.<fnref n="4" />

We have thus argued that, under the ideal scenario and in the weak dissipation regime, the exponential growth of the K-complexity saturates at $\xi$:

$$ K(t) \sim \begin{cases} e^{2\alpha t} & t \lesssim t_*, \\ \xi & t \gtrsim t_* \end{cases}, \quad t_* \sim \frac{1}{2\alpha} \ln \left( \frac{2\alpha}{\chi\mu} \right). \qquad (6.7) $$

Here we also identified the saturation time scale, which diverges as a log of the inverse dissipation strength as $\mu \to 0$ — a behavior reminiscent of the scrambling time with respect to a finite system. At strong dissipation, $t_* \to 1$ and $\xi \to 1$, the exponential growth regime of $K(t)$ disappears.

It is interesting to note that, the “localization” of the operator wavefunction in the region $n \lesssim \xi \propto \alpha/\mu$ on the Krylov chain protects itself against any deviation from the ideal scenario, provided that happens only at $n \gg \xi$; we have seen in section 5.2 that this is indeed the case for the dissipative SYK at finite $q$. It is reasonable to expect this to be general, as the deviation from the ideal scenario necessitates dissipation. Therefore, we believe the K-complexity saturation to be a robust phenomenon in similar situations.

## 6.2 Exact solution

The above argument, although heuristic, is applicable to any $a_n$ and $b_n$ with the asymptotics (6.1). Here, we check the conclusion with an exact solvable set of coefficients:

$$ b_n^2 = (1 - u^2)n(n - 1 + \eta), \quad a_n = iu(2n + \eta). \qquad (6.8) $$

This is precisely the recursion coefficients of the Meixner polynomials used in ref. [1] (appendix D, eq. (D10), upon identifying $\delta = iu$; in that work, only $\delta = 0$ was useful). These Lanczos coefficients have the asymptotic behavior (6.1), with

$$ \alpha = 1 - u^2, \quad \chi\mu = 2u. \qquad (6.9) $$

Varying $u \in (0, 1)$, we have obtained any relative dissipation strength, $\xi^{-1} = \chi\mu/2\alpha \in (0, \infty)$ (meanwhile $\eta > 0$ does not affect the asymptotic growth rates, but only the offsets). Following closely [1, eq. (D11) and (D9)], it is straightforward to show that

$$ \varphi_n(t) = \frac{\operatorname{sech}(t)^\eta}{(1 + u \tanh(t))^\eta} \times (1 - u^2)^{\frac{n}{2}} \sqrt{\frac{(\eta)_n}{n!}} \left( \frac{\tanh(t)}{1 + u \tanh(t)} \right)^n. \qquad (6.10) $$

Here $(\eta)_n = \eta(\eta+1)\dots(\eta+n-1)$ is the Pochhammer symbol. Since we will use the normalized wavefunction to define the K-complexity, the first line of (6.10) (an $n$ independent pre-factor) can be ignored. In figure 4(a), we plotted a few wavefunction snapshots from the above exact solution for $\eta = 1.5$. We see qualitatively that time evolution is indeed described by the above general argument. Curiously, the coefficients (6.8) have the same

<footnote n="4">⁴To be precise, this is true modulo the global normalization, which decays exponentially in the stationary regime; the decay rate depends on details at $n$ small. But at $n$ large, the spatial profile is still given by (6.6).</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 4. (a) Snapshots of the wave-function on the Krylov chain, from the exact solution (6.8), (6.10) with $\eta = 1.5$ and $u = 0.01$. For display, we rescaled the wavefunctions by a $t$-dependent constant; we also interpolated $\varphi_n$ to non-integer values of $n$. At $t \to \infty$, a stationary profile is reached with an exponential tail $\propto e^{-n/\xi}$ with $\xi = 1/u$, indicated by the dashed curve. (b) Growth of the K-complexity for various dissipation strength ($u$), from the exact solution (6.15). The exponential growth, which lasts indefinitely in the closed system ($u = 0$), saturates at $t \sim \ln(1/u)$, or $K(t) \sim 1/u$, in the presence of dissipation.</figcaption></figure>

structure of SL(2, $\mathbb{R}$), after identifying the coefficients [30] and the highest weight state. Hence, the amplitude (6.10) is fixed by the SL(2, $\mathbb{R}$) symmetry, which possibly controls the complexity growth [31]. As we see, in our case, the SL(2, $\mathbb{R}$) structure is preserved, yet the complexity is suppressed. The suppression entirely comes due to the imaginary part of the coefficients $a_n$.<fnref n="5" />

We can make further quantitative comparisons, focusing on the tail of the wave function. It is not hard to see that at $t \to \infty$, the wavefunction has the following stationary $n$-dependence for $u > 0$:

$$ \varphi_n(t \to \infty) \approx \left( \frac{\sqrt{1-u^2}}{1+u} \right)^n n^{\frac{\eta-1}{2}}. \qquad (6.11) $$

We find an exponential decay.<fnref n="6" /> At small $u$ (weak dissipation), its inverse width

$$ \ln \left( \frac{1+u}{\sqrt{1-u^2}} \right) = u + O(u^3), $$

agrees with the prediction of (6.6) applied to (6.9), which gives also $\xi^{-1} = \chi\mu/(2\alpha) = u + O(u^3)$. A similar analysis at finite $t$ gives

$$ \varphi_n(t) \sim e^{-n/\xi(u,t)} n^{\frac{\eta-1}{2}}, \qquad (6.12) $$

where

$$ \xi(u,t)^{-1} = u + 2e^{-2t} + O(e^{-4t}, e^{-2t}u, u^2). \qquad (6.13) $$

<footnote n="5">⁵We thank Pawel Caputa for pointing this out.</footnote>

<footnote n="6">⁶The power-law correction comes from $\sqrt{(\eta)_n/n!}$, after using the asymptotics $\Gamma(z+p)/\Gamma(z+q) \sim z^{p-q}$ as $z \to \infty$.</footnote>
The saturation time is given by equating the two terms; we thus find $t = \ln(1/u)/2 + O(1)$ which agrees with (6.7) at small $u$.

The exact solution above allows us to compute the K-complexity analytically. Indeed, denoting

$$e^{y_0} = (1 - u^2) \left( \frac{\tanh(t)}{1 + u \tanh(t)} \right)^2, \qquad (6.14)$$

it is not hard to see from (6.10) that

$$ \begin{aligned} K(t) &= \frac{1}{\mathcal{Z}} \sum_n |\varphi_n(t)|^2 = \partial_y \ln \left( \sum_n e^{y n \frac{(\eta)_n}{n!}} \right) \Big|_{y=y_0} \\ &= \partial_y \ln \left( (1 - e^y)^{-\eta} \right) \Big|_{y=y_0} \\ &= \frac{\eta (1 - u^2) \tanh^2(t)}{1 + 2u \tanh(t) - (1 - 2u^2) \tanh^2(t)}. \end{aligned} \qquad (6.15) $$

See figure 4(b) for a plot. Expanding this around $u = 0$ (weak dissipation) at fixed $t$ gives

$$K(t) = \eta \left[ \sinh^2(t) - 2u \sinh^3(t) \cosh(t) + O(u^2) \right], \qquad (6.16)$$

where we recover the exponential growth $K(t) \sim e^{2\alpha t}$ (with $\alpha = 1$) in the closed system, as well as the first correction due to dissipation. Now, the latter will dominate at $t \sim t_* \sim \ln(1/u)/2$, rendering this expansion useless. Meanwhile, expanding around $t \to \infty$ at fixed $u > 0$, we see that the K-complexity goes to a constant:

$$K(t \to \infty) = \frac{\eta}{2u} - \frac{\eta}{2}. \qquad (6.17)$$

It is proportional to $\xi \sim 1/u$ at small $u$, and the proportionality constant is controlled by $\eta$. Hence, the exact calculation confirms the conclusions (6.7) of the general argument.

Similarly, we can compute the normalized variance [32, 33] (note that we define it again using the normalized wavefunction):

$$ \begin{aligned} (n^2)_c := \frac{1}{\mathcal{Z}} \sum_n |\varphi_n(t)|^2 (n - K(t))^2 &= \partial_y^2 \ln \left( \sum_n e^{y n \frac{(\eta)_n}{n!}} \right) \Big|_{y=y_0} \\ &= \frac{\eta (1 - u^2) \tanh^2(t) (u \tanh(t) + 1)^2}{(1 + 2u \tanh(t) - (1 - 2u^2) \tanh^2(t))^2}. \end{aligned} \qquad (6.18) $$

It behaves as $\eta e^{4t}$ in the growth regime and saturates at $\eta/(4u^2)$ in the late time regime. In both cases, the position fluctuation $\sqrt{(n^2)_c}$ is comparable to its average $K(t)$.

It is interesting to remark that the Lanczos coefficient (4.2), (4.3) of the large $q$ SYK corresponds exactly — up to subleading terms in a $1/q$ expansion — to the exact solvable ones (4.3) upon identifying

$$\eta = 2/q, \quad 2u = \tilde{\mu}, \quad \mathcal{J}^2 = 1 - u^2. \qquad (6.19)$$

As such, our results above apply to the large $q$ SYK; in particular, we may recover (4.20) and (4.31). The K-complexity in this model is $\propto 1/q$, as most of the wave function is localized at the origin.
## 6.3 Bound on chaos

We now come back to discussing a general consequence of the saturated K-complexity growth (6.7). Now, under the ideal scenario, the K-complexity is computed with respect to the closed-system Krylov basis (which is also approximately the basis generated by the Arnoldi iteration). Thus, the argument leading to the bound on chaos in ref. [1] still applies and entails that the exponential growth of the normalized out-of-time-order correlator (OTOC) must saturate. More precisely and concretely, in the dissipative SYK model, we expect (and have proven at large $q$) that there exists a $t$-independent constant such that

$$ \frac{1}{\mathcal{Z}} \sum_{j=1}^{N} \mathrm{Tr}[\rho_{\infty}\{\psi_j, \psi_1(t)\}\{\psi_j, \psi_1(t)\}^{\dagger}] \leq CK(t), \quad (6.20) $$

where $K(t)$ is given by (6.7), with $\chi$ and $\alpha$ calculable from the model parameters, and $\mathcal{Z} = \mathrm{Tr}[\rho_{\infty}\psi_1(t)^{\dagger}\psi_1(t)]$ is the norm of the evolved operator. In particular, in the presence of any dissipation, however small it is, the OTOC at the left-hand side cannot grow indefinitely even in the $N \to \infty$ limit. In particular, there is no nonzero critical value of $\mu_c$ below which the OTOC is in a “chaotic” phase, in contrast to the models studied recently [12, 17]. Indeed, the dissipation term in our case has a more severe effect on the terms in $\psi_1(t)$ with large size (long Majorana strings), which dominate the OTOC. They are simply eliminated (with a rate proportional to the length); meanwhile, in the models of [12, 17], the dissipation reduces their sizes without eliminating them altogether.

## 6.4 Holographic interpretation

Finally, we speculate some heuristic explanations of our results from the dual gravity picture. As argued in [50], the SYK model with the jump operators given by (3.3) is dual to a Keldysh wormhole. The gravitational geometry of such a wormhole is far from clear. The configuration is valid in real-time and is fundamentally different from the traversable [51] or Euclidean wormhole [59]. However, we can still take inspiration from the bulk picture of the growth of operator size [60, 61]. The scrambling time in the dissipative SYK appears to be smaller than in the non-dissipative SYK. Increasing the coupling decreases the scrambling time, encountering a possible competition between the scrambling and decoherence [7]. In the dual picture, we might think that the particle falling on one side can be reached on the other side due to the coupling like eq. (3.6). However, there is still a finite decay rate even at zero dissipation [48, 50], which is possibly related to the length of the wormhole. The dissipation effectively increases the coupling between the two sides of the wormhole, thus signaling a quicker scrambling. It will be interesting to understand this in a much broader sense, especially computing the complexity in the Jackiw-Teitelboim (JT) gravity [22]. We leave a detailed discussion in future studies.
# 7 Conclusion

In this paper, we aim to extend the universal operator growth hypothesis in open quantum systems. Our model of study is the dissipative $q$-body SYK model where the interactions with the environment are modeled by the Lindbladian construction. One of the particular motivations is the large $q$ limit, where the quantities of interest can be computed analytically. In particular, the “operator size concentration” allows a simple proof to extract two sets of Lanczos coefficients, which can also be exactly computed from the moment method. Out of two sets, only one of the sets is sensitive to dissipation while the other keeps track of the chaotic nature of the Hamiltonian. The asymptotic growth of both sets of coefficients is linear, which can be considered a generalized operator growth hypothesis in the dissipative systems. This result is consistent with the semi-analytics in $q$ by solving the Schwinger-Dyson equations and explicit numerical implementation of Arnoldi iteration. The dissipative nature of the environment suppresses the complexity compared to its exponential growth in a closed system. This result corroborates the finds of SYK for finite $N$ (thus finite $q$) and spinless fermionic systems [20]. We believe that the exponential growth and the consecutive saturation at “scrambling time” is a universal feature to hold any dissipative systems. Finally, we argue some holographic interpretation of the complexity growth of the dissipative SYK, known to be dual to the Keldysh wormhole [50].

While the asymptotic growth of $b_n$s has been known for integrable systems, it is an open question to ask about the asymptotic growth of the $a_n$s in generic integrable systems. Some relevant studies have been performed in [19], and the growth seems to be not affected due to the integrability of the system. We believe a careful study is required. It will also be worth implementing the numerical bi-Lanczos algorithm [62], where the “bra” and “ket” vectors in the doubled Hilbert space will evolve differently. The Lindbladian in this biorthonormal basis is expected to reproduce the analytical results in the large $N$ and large $q$ limit. Finally, it is also exciting to see the growth of the complexity with random jump operators [49], and especially for the dissipative dynamics of complex fermions, where the latter might be equivalent to two non-Hermitian SYK [63, 64]. This might provide a general perspective to establish a “universal” conclusion of our results. A larger and wider motivation is to understand the chaos bound [45] in dissipative systems, provided it is well-defined in such a scenario. We hope to return to some of these questions in future works.

# Acknowledgments

X.C. thanks Thomas Scaffidi and Zihao Qi for useful discussions, in particular for pointing out ref. [20]. P.N. thanks Pawel Caputa, Tokiro Numasawa, and Shinsei Ryu for the fruitful discussions. P.N. would like to thank the organizers of Kavli Asian Winter School (KAWS) 2023 at Daejeon, Korea, and the young researchers' workshop of Extreme Universe (ExU) collaboration at Nagoya University, Japan, where the results were presented. B.B. is supported by the Ministry of Human Resource Development, Government of India through the Prime Ministers' Research Fellowship. The work of P.N. is supported by the JSPS Grant-in-Aid for Transformative Research Areas (A) “Extreme Universe” No. 21H05190.
## A Derivation of eq. (3.4)

In this appendix, we outline a derivation of eq. (3.4). Notice that our jump operators (3.3) are always fermionic. Depending on the string operators, we have the following two cases:

**Case 1:** when the string operators are fermionic, we have to use the negative sign in the Lindbladian (2.4). In this case, the string of operators has to be even to be a fermionic operator. With the jump operators (3.3), we have

$$ \mathcal{L}_D \psi_{i_1} \cdots \psi_{i_s} = i\mu \sum_{k=1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k + \frac{i\mu}{2} \sum_{k=1}^{N} \left\{ \frac{I}{2}, \psi_{i_1} \cdots \psi_{i_s} \right\}, \quad (A.1) $$

for odd $s$. Here we have used $\psi_k^\dagger = \psi_k$ and $\psi_k^2 = I/2$, where $I$ is the identity matrix. The second summation is easier to compute and it gives

$$ \frac{i\mu}{2} \sum_{k=1}^{N} \left\{ \frac{I}{2}, \psi_{i_1} \cdots \psi_{i_s} \right\} = \frac{i\mu N}{2} \psi_{i_1} \cdots \psi_{i_s}. \quad (A.2) $$

The first summation is more involved. To simply it, we break it into two parts

$$ i\mu \sum_{k=1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k = i\mu \sum_{k=1}^{i_s} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k + i\mu \sum_{k=i_s+1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k. $$

Now the first sum is up to the string length $s$. It requires $(s-1)$ shifts to rearrange according to the string length. On the other hand, the second sum involves $(N-s)$ terms and requires $s$ shifts. Hence, the result is

$$ \begin{aligned} i\mu \sum_{k=1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k &= \frac{(-1)^{s-1}}{2} i\mu s \psi_{i_1} \cdots \psi_{i_s} + \frac{(-1)^s}{2} i\mu(N-s) \psi_{i_1} \cdots \psi_{i_s}, \\ &= \frac{(-1)^s}{2} i\mu(N-2s) \psi_{i_1} \cdots \psi_{i_s}. \end{aligned} $$

Since $s$ is odd, we readily have

$$ i\mu \sum_{k=1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k = -\frac{1}{2}i\mu(N-2s)\psi_{i_1} \cdots \psi_{i_s}. $$

Combining with (A.2), we immediately get (3.4). Notice that, as $N$ is even, the r.h.s. of the above equation cannot be zero, it can be positive or negative depending on the number of fermions and the string length.

**Case 2:** when the string operator is bosonic, we have to use the plus sign in the Lindbladian (2.4). In this case, the string of operators has to be even to be a bosonic operator. With the same jump operators (3.3), now we have

$$ \mathcal{L}_D \psi_{i_1} \cdots \psi_{i_s} = -i\mu \sum_{k=1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s})\psi_k + \frac{i\mu}{2} \sum_{k=1}^{N} \left\{ \frac{I}{2}, \psi_{i_1} \cdots \psi_{i_s} \right\}, \quad (A.3) $$
for even $s$. The previous computations will exactly follow. However, now we have

$$i\mu \sum_{k=1}^{N} \psi_k(\psi_{i_1} \cdots \psi_{i_s}) \psi_k = \frac{1}{2} i\mu(N-2s) \psi_{i_1} \cdots \psi_{i_s}. \quad (A.4)$$

Moreover, the eq. (A.2) will still hold for even $s$. Hence, combining (A.2) with (A.3), we find (3.4). Notice that, when $N = 2s$, the first summation in (A.3) will vanish as evident from eq. (A.4). This can be straightforwardly checked by choosing a particular $N$ and $s$.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. SCOAP³ supports the goals of the International Year of Basic Sciences for Sustainable Development.

## References

[1] D.E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi and E. Altman, *A Universal Operator Growth Hypothesis*, Phys. Rev. X 9 (2019) 041017 [arXiv:1812.08657] [INSPIRE].

[2] H.-P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems*, Oxford University Press, Oxford, U.K. (2007).

[3] G. Lindblad, *On the Generators of Quantum Dynamical Semigroups*, Commun. Math. Phys. 48 (1976) 119 [INSPIRE].

[4] V. Gorini, A. Kossakowski and E.C.G. Sudarshan, *Completely Positive Dynamical Semigroups of N Level Systems*, J. Math. Phys. 17 (1976) 821 [INSPIRE].

[5] R. Grobe, F. Haake and H.-J. Sommers, *Quantum Distinction of Regular and Chaotic Dissipative Motion*, Phys. Rev. Lett. 61 (1988) 1899 [INSPIRE].

[6] G. Akemann, M. Kieburg, A. Mielke and T. Prosen, *Universal Signature from Integrability to Chaos in Dissipative Open Quantum Systems*, arXiv:1910.03520 [DOI:10.1103/PhysRevLett.123.254101].

[7] Z. Xu, A. Chenu, T. Prosen and A. del Campo, *Thermofield dynamics: Quantum Chaos versus Decoherence*, Phys. Rev. B 103 (2021) 064309 [arXiv:2008.06444] [INSPIRE].

[8] J. Li, T. Prosen and A. Chan, *Spectral Statistics of Non-Hermitian Matrices and Dissipative Quantum Chaos*, Phys. Rev. Lett. 127 (2021) 170602 [arXiv:2103.05001] [INSPIRE].

[9] A.M. García-Garcí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 (2022) 021040 [arXiv:2110.03444] [INSPIRE].

[10] A.S. Matsoukas-Roubeas, F. Roccati, J. Cornelius, Z. Xu, A. Chenu and A. del Campo, *Non-Hermitian Hamiltonian deformations in quantum mechanics*, JHEP 01 (2023) 060 [arXiv:2211.05437] [INSPIRE].

[11] K. Kawabata, A. Kulkarni, J. Li, T. Numasawa and S. Ryu, *Symmetry of open quantum systems: Classification of dissipative quantum chaos*, arXiv:2212.00605 [INSPIRE].

[12] P. Zhang and Z. Yu, *Dynamical Transition of Operator Size Growth in Open Quantum Systems*, arXiv:2211.03535 [INSPIRE].
[13] S. Omanakuttan, K. Chinni, P.D. Blocher and P.M. Poggi, *Scrambling and quantum chaos indicators from long-time properties of operator distributions*, arXiv:2211.15872 [INSPIRE].

[14] P. Zhang and Y. Gu, *Operator Size Distribution in Large N Quantum Mechanics of Majorana Fermions*, arXiv:2212.04358 [INSPIRE].

[15] P. Zanardi and N. Anand, *Information scrambling and chaos in open quantum systems*, Phys. Rev. A 103 (2021) 062214 [arXiv:2012.13172] [INSPIRE].

[16] T. Schuster and N.Y. Yao, *Operator Growth in Open Quantum Systems*, arXiv:2208.12272 [INSPIRE].

[17] Z. Weinstein, S.P. Kelly, J. Marino and E. Altman, *Scrambling Transition in a Radiative Random Unitary Circuit*, arXiv:2210.14242 [INSPIRE].

[18] S.V. Syzranov, A.V. Gorshkov and V. Galitski, *Out-of-time-order correlators in finite open systems*, Phys. Rev. B 97 (2018) 161114 [arXiv:1704.08442] [INSPIRE].

[19] A. Bhattacharya, P. Nandy, P.P. Nath and H. Sahu, *Operator growth and Krylov construction in dissipative open quantum systems*, JHEP 12 (2022) 081 [arXiv:2207.05347] [INSPIRE].

[20] C. Liu, H. Tang and H. Zhai, *Krylov Complexity in Open Quantum Systems*, arXiv:2207.13603 [INSPIRE].

[21] E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, *Operator complexity: a journey to the edge of Krylov space*, JHEP 06 (2021) 062 [arXiv:2009.01862] [INSPIRE].

[22] S.-K. Jian, B. Swingle and Z.-Y. Xian, *Complexity growth of operators in the SYK model and in JT gravity*, JHEP 03 (2021) 014 [arXiv:2008.12274] [INSPIRE].

[23] J.L.F. Barbón, E. Rabinovici, R. Shir and R. Sinha, *On The Evolution Of Operator Complexity Beyond Scrambling*, JHEP 10 (2019) 264 [arXiv:1907.05393] [INSPIRE].

[24] A. Dymarsky and A. Gorsky, *Quantum chaos as delocalization in Krylov space*, Phys. Rev. B 102 (2020) 085137 [arXiv:1912.12227] [INSPIRE].

[25] X. Cao, *A statistical mechanism for operator growth*, J. Phys. A 54 (2021) 144001 [arXiv:2012.06544] [INSPIRE].

[26] A. Dymarsky and M. Smolkin, *Krylov complexity in conformal field theory*, Phys. Rev. D 104 (2021) L081702 [arXiv:2104.09514] [INSPIRE].

[27] P. Caputa, J.M. Magan and D. Patramanis, *Geometry of Krylov complexity*, Phys. Rev. Res. 4 (2022) 013041 [arXiv:2109.03824] [INSPIRE].

[28] E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, *Krylov localization and suppression of complexity*, JHEP 03 (2022) 211 [arXiv:2112.12128] [INSPIRE].

[29] B. Bhattacharjee, X. Cao, P. Nandy and T. Pathak, *Krylov complexity in saddle-dominated scrambling*, JHEP 05 (2022) 174 [arXiv:2203.03534] [INSPIRE].

[30] V. Balasubramanian, P. Caputa, J.M. Magan and Q. Wu, *Quantum chaos and the complexity of spread of states*, Phys. Rev. D 106 (2022) 046007 [arXiv:2202.06957] [INSPIRE].

[31] N. Hörnedal, N. Carabba, A.S. Matsoukas-Roubeas and A. del Campo, *Ultimate Speed Limits to the Growth of Operator Complexity*, Commun. Phys. 5 (2022) 207 [arXiv:2202.05006] [INSPIRE].

[32] P. Caputa and S. Datta, *Operator growth in 2d CFT*, JHEP 12 (2021) 188 [Erratum ibid. 09 (2021) 113] [arXiv:2110.10519] [INSPIRE].
[33] B. Bhattacharjee, P. Nandy and T. Pathak, *Krylov complexity in large-q and double-scaled SYK model*, arXiv:2210.02474 [YITP-22-106] [INSPIRE].

[34] B. Bhattacharjee, S. Sur and P. Nandy, *Probing quantum scars and weak ergodicity breaking through quantum complexity*, Phys. Rev. B 106 (2022) 205150 [arXiv:2208.05503] [INSPIRE].

[35] V. Balasubramanian, J.M. Magan and Q. Wu, *A Tale of Two Hungarians: Tridiagonalizing Random Matrices*, arXiv:2208.08452 [INSPIRE].

[36] E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, *Krylov complexity from integrability to chaos*, JHEP 07 (2022) 151 [arXiv:2207.07701] [INSPIRE].

[37] S. He, P.H.C. Lau, Z.-Y. Xian and L. Zhao, *Quantum chaos, scrambling and operator growth in $T\bar{T}$ deformed SYK models*, JHEP 12 (2022) 070 [arXiv:2209.14936] [INSPIRE].

[38] S. Guo, *Operator growth in $SU(2)$ Yang-Mills theory*, arXiv:2208.13362 [INSPIRE].

[39] S. Sachdev and J. Ye, *Gapless spin fluid ground state in a random, quantum Heisenberg magnet*, Phys. Rev. Lett. 70 (1993) 3339 [cond-mat/9212030] [INSPIRE].

[40] J. Maldacena and D. Stanford, *Remarks on the Sachdev-Ye-Kitaev model*, Phys. Rev. D 94 (2016) 106002 [arXiv:1604.07818] [INSPIRE].

[41] J. Maldacena, D. Stanford and Z. Yang, *Conformal symmetry and its breaking in two dimensional Nearly Anti-de-Sitter space*, Prog. Theor. Exp. Phys. 2016 (2016) 12C104 [arXiv:1606.01857] [INSPIRE].

[42] A. Kitaev, *A simple model of quantum holography. Part 1*, talk given at KITP, Santa Barbara, CA, U.S.A., 6 April–2 July 2015 and online at https://online.kitp.ucsb.edu/online/joint98/kitaev/.

[43] A. Kitaev, *A simple model of quantum holography. Part 2*, talk given at KITP, Santa Barbara, CA, U.S.A., 6 April–2 July 2015 and online at https://online.kitp.ucsb.edu/online/entangled15/kitaev2/.

[44] J.S. Cotler et al., *Black Holes and Random Matrices*, JHEP 05 (2017) 118 [Erratum ibid. 09 (2018) 002] [arXiv:1611.04650] [INSPIRE].

[45] J. Maldacena, S.H. Shenker and D. Stanford, *A bound on chaos*, JHEP 08 (2016) 106 [arXiv:1503.01409] [INSPIRE].

[46] C. Lanczos, *An iteration method for the solution of the eigenvalue problem of linear differential and integral operators*, J. Res. Natl. Bur. Stand. B 45 (1950) 255 [INSPIRE].

[47] W.E. Arnoldi, *The principle of minimized iterations in the solution of the matrix eigenvalue problem*, Quart. Appl. Math. 9 (1951) 17.

[48] A. Kulkarni, T. Numasawa and S. Ryu, *Lindbladian dynamics of the Sachdev-Ye-Kitaev model*, Phys. Rev. B 106 (2022) 075138 [arXiv:2112.13489] [INSPIRE].

[49] L. Sá, P. Ribeiro and T. Prosen, *Lindbladian dissipation of strongly-correlated quantum matter*, Phys. Rev. Res. 4 (2022) L022068 [arXiv:2112.12109] [INSPIRE].

[50] A.M. García-García, L. Sá, J.J.M. Verbaarschot and J.P. Zheng, *Keldysh Wormholes and Anomalous Relaxation in the Dissipative Sachdev-Ye-Kitaev Model*, arXiv:2210.01695 [INSPIRE].

[51] J. Maldacena and X.-L. Qi, *Eternal traversable wormhole*, arXiv:1804.00491 [INSPIRE].
[52] V. Viswanath and G. Müller, *The Recursion Method: Application to Many Body Dynamics*, in *Lecture Notes in Physics Monographs*, Springer (1994).

[53] K. Kawabata, A. Kulkarni, J. Li, T. Numasawa and S. Ryu, *Dynamical quantum phase transitions in SYK Lindbladians*, arXiv:2210.04093 [insPIRE].

[54] M.-D. Choi, *Completely positive linear maps on complex matrices*, Linear Algebra Appl. **10** (1975) 285.

[55] A. Jamiołkowski, *Linear transformations which preserve trace and positive semidefiniteness of operators*, Rep. Math. Phys. **3** (1972) 275.

[56] D.A. Roberts, D. Stanford and A. Streicher, *Operator growth in the SYK model*, JHEP **06** (2018) 122 [arXiv:1802.02633] [insPIRE].

[57] D. Knuth, *The online encyclopedia of integer sequences: A101280*, (2005) https://oeis.org/A101280.

[58] X.-L. Qi and A. Streicher, *Quantum Epidemiology: Operator Growth, Thermal Effects, and SYK*, JHEP **08** (2019) 012 [arXiv:1810.11958] [insPIRE].

[59] A.M. García-García and V. Godet, *Euclidean wormhole in the Sachdev-Ye-Kitaev model*, Phys. Rev. D **103** (2021) 046014 [arXiv:2010.11633] [insPIRE].

[60] L. Susskind, *Why do Things Fall?*, arXiv:1802.01198 [insPIRE].

[61] L. Susskind and Y. Zhao, *Complexity and Momentum*, JHEP **03** (2021) 239 [arXiv:2006.03019] [insPIRE].

[62] M. Grüning, A. Marini and X. Gonze, *Implementation and testing of Lanczos-based algorithms for Random-Phase Approximation eigenproblems*, Comput. Mater. Sci. **50** (2011) 2148 [arXiv:1102.3909].

[63] A.M. García-Garcí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*, arXiv:2211.01650 [insPIRE].

[64] W. Cai, S. Cao, X.-H. Ge, M. Matsumoto and S.-J. Sin, *Non-Hermitian quantum system generated from two coupled Sachdev-Ye-Kitaev models*, Phys. Rev. D **106** (2022) 106010 [arXiv:2208.10800] [insPIRE].