Operator dynamics in Lindbladian SYK: a Krylov complexity perspective

We use Krylov complexity to study operator growth in the $q$-body dissipative SYK model, where the dissipation is modeled by linear and random $p$-body Lindblad operators. In the large $q$ limit, we analytically establish the linear growth of two sets of coefficients for any generic jump operators. We numerically verify this by implementing the bi-Lanczos algorithm, which transforms the Lindbladian into a pure tridiagonal form. We find that the Krylov complexity saturates inversely with the dissipation strength, while the dissipative timescale grows logarithmically. This is akin to the behavior of other $\mathfrak{q}$-complexity measures, namely out-of-time-order correlator (OTOC) and operator size, which we also demonstrate. We connect these observations to continuous quantum measurement processes. We further investigate the pole structure of a generic auto-correlation and the high-frequency behavior of the spectral function in the presence of dissipation, thereby revealing a general principle for operator growth in dissipative quantum chaotic systems.

local operator grows under the time evolution by the Hamiltonian of the system.Several ways of characterizing such growth have been proposed in recent years, namely the operator size distribution [1,2], out-of-time-ordered correlator (OTOC) [3], and Krylov complexity [4].These approaches are indirect since they require a probe operator and contrast with more direct probes such as level statistics [5,6] and spectral form factor (SFF) [7].The latter is particularly useful in the semiclassical theory of gravity, where a dip, ramp, and plateau are observed and feature behavior similar to underlying random matrix universality at late times.On the other hand operator growth turns out to be useful in probing the black hole interior.Particularly the increasing momentum of a particle in a near AdS 2 black hole is reflected as an operator growth in the dual geometry [8], often studied under the shadow of operator complexity [9,10].
In this paper, our interest is in studying Krylov complexity, which measures operator growth on a special basis known as the Krylov basis.The basis is formed by an iterative Gram-Schmidt-like orthonormalization procedure known as the Lanczos algorithm [11,12].An output of this algorithm is the Lanczos coefficients, which usually show distinctive features for integrable and non-integrable systems, albeit in some special cases [13,14].The universal operator growth hypothesis, proposed in [4] states that the Lanczos coefficients can at most grow linearly in their indices, and chaotic systems exhibit this fastest linear growth.This linear growth is consistent with the chaos bound [3] and has attracted substantial studies in recent times , along with its cousins, circuit complexity and holographic complexity [39].
On the other hand, open quantum systems are quantum systems that interact with their environment, which causes decoherence and dissipation.These effects are common in nature and have practical implications for quantum many-body systems.In particular, the study is important in the context of the black hole information problem in AdS/CFT correspondence, where the Hawking radiation is collected in a bath that is attached to an AdS black hole [40].In recent times, this has led to a surge of studies on the operator dynamics in open quantum systems, from quantum many-body systems [41][42][43][44] to quantum field theory and holography [45].
However, the operator evolution in an open system is drastically different compared to a closed system.The primary reason is the interaction with the environment, which makes the whole evolution non-unitary.This poses a challenge for studying operator growth in open systems, as the usual Lanczos algorithm does not work [21].The problem can be circumvented by applying more generic algorithms such as Arnoldi iteration [46] or the bi-Lanczos algorithm.A generic study of Krylov complexity was initiated in [21][22][23] using such algorithms.Especially, the authors of this paper initiated a study in the dissipative Sachdev-Ye-Kitaev (SYK) model [47] using Arnoldi iteration and motivated some universal aspects of the growth and saturation of Krylov complexity [22].This study is relevant since this particular dissipative SYK model can be interpreted as two non-Hermitian SYK models coupled by a Keldysh wormhole [48], suggesting that generic open quantum systems may have a holographic dual that involves a wormhole or a similar structure.However, a special form of dissipation is chosen, which is simple and thus not generic.Therefore, it is still an open question as to how robust the results are for other forms of dissipation.
In this paper, we partially answer this question by studying a large class of dissipators.By choosing a p-body Lindblad dissipator with Gaussian strength, in the large N and the large q limit, we analytically show that both the diagonal and off-diagonal coefficients of the Lindbladian matrix exhibit asymptotically linear growth, consistent with the observation made in spin chains [23].This is further supported by the results of the bi-Lanczos algorithm in the finite N and finite q SYK in the appropriate regime, modulo the finite-size effects.The resulting logarithmic timescale of dissipation and the saturation of Krylov complexity are found to be fairly general and independent of the choice of the form of the dissipation.We find that this growth and saturation are also reflected in the behavior of OTOC and operator size, which is supposed to construct a larger class of q-complexity measures [4].The saturation can also be interpreted as a result of continuous measurement by the environment itself.Finally, we also provide a generic notion of the pole structure of auto-correlation and the high-frequency behavior of the spectral function in the presence of dissipation.This leads us to motivate an operator growth hypothesis in generic open quantum systems.
Our paper is structured as follows.In section 2, we introduce our system and the generic form of the dissipation.Section 3 introduces the machinery of the bi-Lanczos algorithm and the general version of the Krylov complexity in dissipative systems.In section 4, we provide an analytical derivation of the linear growth of the diagonal coefficients of the Lindbladian matrix for any generic dissipation which we numerically confirm by implementing the bi-Lanczos algorithm in section 5. Based on the above results, section 6 motivates some universal aspects of Krylov complexity and its associated quantities.Finally, in section 7, we derive the expression of OTOC for a 1-body and general p-body fermionic initial operators.We conclude in section 8 with a brief summary and future outlook.Appendices consist of some further results and detailed derivation which we omit in the main text.
2 System and the environment: The Lindbladian The prototypical example of an open system consists of a system which is interacting with a dissipative environment.Our system under study is the Sachdev-Ye-Kitaev (SYK) model [49,50].This model has garnered much attention in recent times, especially being maximally chaotic [7,51], and sharing the same Schwarzian action as Jackiw-Teitelboim (JT) gravity in low temperatures which elevates it as a toy model for holography.In addition, the generic model is known to be maximally chaotic, satisfying the Maldacena-Shenker-Stanford (MSS) bound [3].From a condensed matter physics perspective, it provides detailed insights into non-Fermi liquids and strange metals [52].
The q-body Sachdev-Ye-Kitaev (SYK) model consists of N Majorana fermions ψ i , satisfying the Clifford algebra {ψ a , ψ b } = δ ab , where q fermions are interacting at a time.The Hamiltonian is given by [49,50] The random couplings J i 1 ⋯iq are drawn from the Gaussian ensemble with the following mean and variance where J 2 = 21−q qJ 2 .This notation is specifically useful in the large q limit and N → ∞ limit since the model is chaotic and becomes analytically tractable in this limit.As we have stated before, we will treat this as a system and examine various variations of this model in the open system setting in the following subsections.We consider the system to be connected to an environment governed by Markovian dynamics.This regime is defined when the system density matrix ρ(t + dt) = ρ(t) + O(dt) is solely determined by the system density matrix at time t i.e., ρ(t).More broadly, we consider the dissipative mechanism where the information leaks out from the system to the environment such that it never returns to the system at a later time. 1 In other words, our observed time scale of the system dynamics t sys is long compared to the timescale δt E that the environment retains the memory of the information that has been leaked out from the system i.e., t sys ≫ δt E .Under the Born-Markovian approximation, and weak coupling regime, the evolution of the density matrix and any operator can be treated in the realm of Lindbladian formalism [53,54].The density matrix of the system evolves by the master equation where H is the system Hamiltonian.The operators L k are referred to as Lindblad jump operators and they capture the information of the interaction between the system and the environment.In particular, they are made of system operators only and completely lack detailed information about the environment.An arbitrary initial operator O 0 at t = 0 evolves as Here L † o is known as the (adjoint) Lindbladian for the operator2 and acts as Here the "−" sign should be considered in case both L k and O are fermionic [24].One usually express it in a vectorized form after the Choi-Jami lkowski isomorphism [55,56] with the following replacement [57,58] where A and B are any arbitrary operators and vec O is the vectorization of the operator O.This gives the Lindbladian superoperator The notation "≡" indicates that the Lindbladian is represented in a matrix form in the doubled-Hilbert space.In this paper, our system is the SYK Hamiltonian (2.1) and we take the following two classes of jump operators: Class 1: Linear dissipator: We consider an open system version of SYK with the following jump operators [47]: with λ ≥ 0 being the coupling strength between the system and the environment.We often call the full system a dissipative SYK with linear jump operators.This is the simplest version of dissipative SYK, where each of the fermions dissipates at an equal rate.One can solve this model analytically in the large-q limit.This model is particularly useful in being realized as a connected Keldysh wormhole [48].A detailed study of Krylov complexity in this setup was conducted in [22].This analysis can be extended to a dissipation strength of the type V i (instead of √ λ) where V i are Gaussian random complex numbers with zero mean and finite variance.The results are the same under disorder averaging, up to an appropriate identification of the respective parameters.
Class 2: Non-linear dissipator: For this class, we take the p-body jump operators of the following form [59] with the following distribution of V i 1 i 2 ...ip : with V ≥ 0. In other words, the jump operators are p-local and mimic the SYK-like structure.Together with the Hamiltonian (2.1), they dictate the full non-unitary dynamics governed by the Lindbladian (2.7).In particular, the parameter J represents the unitary dynamics while the parameter V breaks it.The p = 1 case without the random average is known as the linear dissipator (class 1 ).The p = 2 case is known as the quadratic dissipator model, previously introduced in [47,60].In the following sections, our interests will be the generic p-body dissipator with a possible emphasis on linear (p = 1) and quadratic (p = 2) dissipator cases particularly. 3efore jumping on to the numerical machinery of the bi-Lanczos algorithm, we briefly discuss a different perspective of the Lindblad (and Lindblad-like) equation [61].Suppose we prepare the system at time t with a density matrix given by ρ(t) and evolve unitarity till time t + δt.The state of the system is ρ(t + δt), which can be written as (2.11) This equation purely comes from the unitary dynamics of the system, namely the Schrodinger equation dρ(t)/dt = −i[H, ρ(t)].However, if we measure the system with a probability P (t + δt) at time t + δt, then the state after the measurement is given by where L k 's are the same quantum jump operators as introduced in (2.5).They act as projector operators satisfying the completeness relation 4 We further expand the probability as where η(t) = δP (t)/δt is the change of measurement probability in unit time and we refer to it as the measurement rate.Plugging (2.11) and (2.13) into the expansion of (2.12), and neglecting O(δt 2 ) terms, we obtain which has a surprisingly similar form to the Lindblad master equation for the density matrix.The anti-commutator part is trivial here due to the completeness relation.
The above expression (2.14) bears a physical significance.In particular, we can directly associate the dissipation strength in (2.8) or (2.10) as the measurement rate η(t) by the environment itself.Depending on the strength of the dissipation, the system can be either in a fully scrambled phase or a purely dissipative phase, opening the possibility of studying the so-called "environment-induced phase transition" [42].

Bi-Lanczos algorithm for open systems
As we briefly discussed in the introduction section, the bi-Lanczos algorithm is a numerical method that can transform a Lindbladian matrix into a tri-diagonal form, which can be used to compute the Lanczos coefficients and the Krylov complexity of an open quantum system.The bi-Lanczos algorithm was first applied to the study of Krylov complexity in [23], where some properties of the coefficients were studied in spin chains.In this section, we will review the bi-Lanczos algorithm and present some more properties of the coefficients, such as their asymptotic behavior.We will then apply the bi-Lanczos algorithm to the dissipative SYK model in later sections, and compare our results with analytical counterparts.

Vectorization and bi-orthonormal vectors
The central idea of the bi-Lanczos algorithm is to construct two sets of bi-orthonormal vectors |p n ⟫ and |q n ⟫ satisfying the following bi-orthonormal condition We use the notation of [23] to denote the bi-Lanczos vectors with "double braces".They are obtained by vectorizing the initial operator.These two bi-orthonormal vectors evolve differently under the Lindblad evolution.In the absence of dissipation, the Lindbladian reduces to the Liouvillian, which is Hermitian and can be recast into a purely tridiagonal form with vanishing diagonal coefficients.The two spaces become conjugate to each other and thus become individually orthonormal [4].Now, we outline the steps of the bi-Lanczos algorithm [23, 63]:5 1. Initialization.
, where O 0 is the initial vector.
For j = 1, 2, . .., we perform the following iterations: (h) Check the convergence and perform the full orthogonalization (FO) procedure, if required.
3. Stop, if c k = 0 for some k.
The algorithm generates three sets of coefficients {a j }, {b j } and {c j }, and two sets of bi-orthogonal vectors {|p j ⟫} and {|q j ⟫}.The full action of this bi-Lanczos basis can be expressed in the following form, which are two sets of three-term recurrences where * denotes the complex conjugate.In other words, we have generated two sets of Krylov spaces, one acts by L o and the other one by L † o : From the recurrence (3.3), it is evident that the procedure of the bi-Lanczos algorithm recasts the Lindbladian into the following tridiagonal form In other words, we have elements in the diagonal and the primary off-diagonal terms only.This is a distinctive feature from the Arnoldi iteration [21,22] which considered only a single set of orthonormal vectors and renders the Lindbladian into an upper-Hessenberg form.However, the methodology to generate the orthonormal vectors is different in the Arnoldi iteration where the individual space is orthonormal.Also, the numerical stability significantly differs in both cases.While the computational cost (time complexity) is higher in the Arnoldi iteration compared to the bi-Lanczos algorithm, the latter might suffer a breakdown i.e., the loss of orthogonality.Such breakdown never occurs in Arnoldi iteration [64].
As we will see in the upcoming sections, our numerical algorithm implies b n = c n .The Lindbladian matrix can be written in tridiagonal form in Krylov (bi-Lanczos) basis in the following form We find that the diagonal coefficients a n are purely imaginary a n = i|a n | and the off-diagonal elements b n are purely real.This implies that the Lindbladian It is generically non-Hermitian.For more details of the properties of these coefficients, see [23] and subsection 3.3.

Operator growth and Krylov complexity
Since the operator evolves with L † o , we expand the time-evolved operator in the bi-Lanczos basis in the following form where φ n are Krylov basis wavefunctions.We slightly changed the notation from [4] to keep the auto-correlation function as φ 0 since our initial vector starts with |p 1 ⟫ (and where φ n ≡ φ n (t) for brevity and φ 0 is the auto-correlation function.For the Lindbladian evolution, it is defined as where O(t) is given by Eq. (2.4), {µ} is the set of the dissipative parameters, and N is the number of degrees of freedom in the system.Eq. (3.10) can be interpreted as a non-Hermitian tight-binding model [24], which we refer to as the non-Hermitian Krylov chain (see Fig. 1).The particle hops from n-th site to (n + 1)-th site with hopping amplitude b n+1 and to (n − 1)-th site with a different amplitude c n , while the amplitude of staying at that particular site is a n+1 .As we will examine later, in all the examples of various versions of dissipative SYK, we numerically find b n = c n for all n.Moreover we find all a n are imaginary, i.e., a n = i|a n |.Hence the Eq.(3.10) simplifies to6 The Krylov complexity is thus defined as the average position of a particle in the non-Hermitian Krylov chain given by where Z = ∑ n |φ n (t)| 2 is the normalization.The probability ∑ n |φ n (t)| 2 < 1 is not conserved due to the unitarity breaking, thus a division is required since the rescaled amplitude φ(t)/ √ Z conserves the probability i.e., ∑ n |φ n (t)/ √ Z| 2 = 1.

Properties of bi-Lanczos coefficients
The bi-Lanczos algorithm generates three sets of coefficients.In general, all the coefficients can be complex numbers.However, as we see the structure and physical properties of a Lindbladian matrix heavily restrict the properties of these coefficients.The density matrix evolution is governed by L o , while the operator evolution is governed by the adjoint L † o [41].This will be reflected transparently in the structure of Lindbladian in the bi-Lanczos basis.Specifically, L † o has positive and purely imaginary diagonal coefficients, while the diagonal coefficients of L o are negative and purely imaginary.In either case, the off-diagonal coefficients are purely real.Together, this makes L † o or L o non-Hermitian.The structure of these elements must be consistent such that the eigenvalues of −iL † o or iL o have to be either purely real positive or complex conjugate in pairs with the real part being positive. 7In all the examples we study, these conditions will be fulfilled.However, before delving into the examples, we state a proposition8 concerning the elements of the bi-Lanczos coefficients.

Proposition 1. The imaginary part of any eigenvalue λ
where K is the Krylov dimension.For a closed system, λ L of L † o is Hermitian, and all a n vanishes, thus (3.14) holds trivially with equality.
Proof: Following [66], we consider a diagonal matrix D, such that we transform the matrix (3.8) into the following form where the sum runs over the corresponding elements (up to the Krylov dimension).
The form of the Lindbladian (3.15) suggests an alternate form of Lanczos coefficients d n ∶= √ b n c n , which was recently advocated in [67].

Lindbladian SYK: analytical treatment
In this section, we provide a simple analytical treatment to show the linear growth of diagonal coefficients.We split the Lindbladian (2.5) into two parts Here we need to choose the "-" sign if both the jump operators and the initial operators are fermionic.Our derivation is based on a property known as the "operator size concentration" [22].It states that the eigenoperators of the dissipative part L † D (ignoring o(1/q) corrections) are given by where c i 1 i 2 ⋯is are some coefficients, and s = n(q − 2) + 1.This property only holds in the large q and large N limit, hence our analytical derivation only holds in this limit.

For linear dissipator
We will present our first example of this model in the open-system setting.This is done by introducing the linear jump operators of the form [47]: with λ ≥ 0. These are the jump operators of class 1, as we have discussed earlier.
For the case of the fermionic jump operator, it was shown in [22] that the dissipative part of the Lindbladian acts linearly i.e., where λ = λq.This immediately gives a n ∼ i λn.Moreover, the advantage of these jump operators is that they allow us to perform an exact analytical calculation in the large-q limit.In particular, one can solve the Schwinger-Dyson equation and obtain the following two-point function [47]: where λ = λq is a redefined coupling in the large-q limit.The parameters α and ℵ are related to the couplings as The closed-system result recovers in the limit λ = 0 and thereby obtaining g(t) = 2 ln sech(J t), which is a known result [68].We can also compute the Lanczos coefficients using the moment method [12].They are given by [22] a n = i λn + o(1/q) , λ ∶= λq , (4.8) Not only does the expression of b n match exactly with the closed-system result [4], but also the expression of a n matches what we obtained from the "operator size concentration" property, with an o(1/q) correction which vanishes in the large q limit.This gives the first hint that the off-diagonal elements of the Lindbladian matrix (3.8) might not depend on the dissipation which is entirely reflected in the diagonal coefficients.However, the linear dissipation (4.3) is very special and it is unclear if the above conclusion is generic for any arbitrary dissipation.Thus, we need to choose a more generic dissipation of class 2 to justify (or falsify) our conclusion.

For random quadratic dissipator
Next, we consider the jump operators which are random and quadratic.This belongs to the class 2 non-linear dissipator with p = 2.In particular, we choose with the following distribution of V ij drawn of random Gaussian ensemble In principle, V can be arbitrary but for our computation, we focus on the weakdissipation regime, which implies J ≫ V . 9This choice is motivated by the fact that we are considering the system dynamics that is Markovian.
Similar to the linear case, we divide the Lindbladian into a Hermitian and a dissipative part.We choose the "+" sign in (4.1) since the jump operators are bosonic.Since our primary concern is the dissipative part L † D , we consider the string ψ i 1 ψ i 2 . . .ψ is with i 1 < i 2 < ⋅ ⋅ ⋅ < i s and attempt to divine its' action on.In fact, the action of the dissipative part of the Lindbladian to a string of length s results in the following proposition: Proposition 2. Under ensemble averaging, the action of the dissipative part of the Lindbladian to a string of length s = n(q − 2) + 1 results in the following expression: with R = M /N .Here ζ ∼ o(1) number and V is given by the ensemble average (4.11).
The proof is given in the Appendix A. The asymptotic linear growth can be deducted In the large q and large N , limit R = M /N becomes the relevant quantity with M being the number of jump operators.We can generalize the above result as a generic p-body dissipator of the form (2.9)-(2.10).However, the computation is tedious and put in the Appendix B. We obtain which is strictly valid in the large N and large q limit.It is easy to see that this reduces to the leading order contribution (4.12) for p = 2.It is straightforward to conclude i.e., the asymptotic growth of the diagonal coefficients is linear and similar to (4.13).

Bi-Lanczos algorithm in Lindbladian SYK
To justify the above analytical results, we resort to the numerical bi-Lanczos algorithm to transform the Lindbladian into a pure tridiagonal form as given by equation (3.8).We separately apply this algorithm for linear, quadratic, and cubic dissipators.

For linear jump operators
Fig. 2 shows the result for N = 20, SYK 4 with 50 Hamiltonian realizations.We can see that the off-diagonal coefficients are unaffected by the dissipation and they   are exactly equal to the closed-system counterparts [4].On the other hand, the dissipation only influences the diagonal coefficients.They are purely imaginary and we can compute the slope of diagonal coefficients as which grows linearly.This slope agrees with the result obtained by Arnoldi iteration in [22], except for some constant shift that depends on the intrinsic nature of the algorithm.Given the set of Lanczos coefficients, Proposition 1 holds which can be checked explicitly.Moreover, as seen from Fig. 3, the slopes of both diagonal and off-diagonal coefficients do not depend on the system size N while their saturation value does.In fact, the saturation linearly increases with the system size N , as shown in the insets of Fig. 3, i.e., for a fixed dissipation strength µ.This finite-size scaling is consistent with previous studies [22,69].However, in the true thermodynamic limit N → ∞, we only observe the asymptotic growth (6.1).

For random quadratic jump operators
With this choice, we perform the bi-Lanczos algorithm for several numbers of Lindblad operators (i.e., different values of M ) with a fixed choice of initial operator O 0 and system size N .The Lanczos coefficients are shown in Fig. 4. We see that the diagonal coefficients |a n | are strongly dependent on the dissipation while the offdiagonal coefficients (b n = c n ) are independent of the dissipation.We also checked that the Proposition 1 remains to hold with the observed set of Lanczos coefficients.We remark on two features of the diagonal coefficients.First, we observe that both the slopes and the saturation values of |a n | increase with M .In Fig. 5a  From our linear dissipator result, we can also understand that increasing the system size N increases the individual saturation value but does not affect the slope.Second, we assume that a n is an asymptotically smooth function of n in the thermodynamic limit.This smoothness behavior is a typical assumption of the operator growth hypothesis for the off-diagonal coefficients [4], although some violations were observed in quantum field theories [29,30] 10 However, in this paper, we continue the smoothness assumption which enables us to define the growth rate of the form with fixed N .In the large N and large q limit, the growth will be asymptotic.In other words, we can write where the proportionality constant c V depends on the dissipation strength V .The second equality comes in a special "double-scaling limit", which is defined as    and N → ∞ keeping R = M /N finite (fixed).Our interest is to find the form of the proportionality constant c V .In principle, it can be either analytically found by computing a two-point function as in (4.5) or numerically by a fitting of the various data of c V .We choose the latter approach.The numerical data obtained by implementing the bi-Lanczos algorithm suggests the following form where κ, β are some real coefficients and can be obtained by fitting the data which is shown in Fig. 5c.Note that this set is obtained for N = 16, and can be improved by increasing N .For our interest, the coefficient κ is irrelevant and we are primarily interested in the exponent β.The fitting suggests β = 2.0046 ≈ 2. Hence, we can write (5.5) of the form This asymptotic growth is consistent with (4.13) and of the form (6.1), while b n follows the same growth of a closed system with coefficient α.

Universal aspects of operator growth in open systems
In this section, we interpret both the analytic and the numerical results into a concise form and discuss some universal aspects of operator growth, generic to any choice of Lindbladian.

An asymptotic growth of Lanczos coefficients and Krylov complexity
The analytical and the numerical analysis motivate us to propose the following sets of asymptotic growth for the Lanczos coefficients in the large-n limit [22] where µ is the generic dissipative parameter, χ is some number which is independent of the dissipation and α captures the information of the Hamiltonian and the initial operator.This form motivates an operator growth hypothesis in open quantum systems.As we have from the previous analysis, the growth of a n is linear, and thus µ ∝ λ for the linear dissipator and µ ∝ RV 2 for the generic p-body dissipator.In either case, the dissipation strength is quadratic due to the simultaneous appearance of L k and L † k in the Lindbladian.One can directly calculate the wavefunctions and the K-complexity from this asymptotic growth (6.1).First, the Krylov basis wavefunctions φ n are given by [22] φ n (t) = sech(γt) η for the exact expression b 2 n = (1−u 2 )γ 2 n(n−1+η) and a n = iuγ(2n+η) which reduces to (6.1) for α 2 = γ 2 (1 − u 2 ) and χµ = 2γu asymptotically for some η ∼ o(1).The Krylov complexity can be straightforwardly computed as [22] The behavior of Krylov complexity is shown in Fig. 6 for different dissipation strengths.
In particular, we observe that the dissipative time scale is logarithmic while the saturation of Krylov complexity scales inversely to the dissipative strength [22]: and finally, reach a value that is independent of the system size.This logarithmic timescale and saturation is also found in operator size distribution [43].
The scaling of the above saturation invites an interpretation from the quantum measurement.Recall (2.14), where the rate of measurement is translated as the dissipation strength in the Markovian approximation.In other words, the jump operators can be interpreted as performing a similar task to measurement operators -the environment makes a continuous measurement through it.However, a significant difference from generic measurement is that here the outcome is unknown to us.However, since the measurement is a non-unitary process, the stronger the measurement rate, the lower the probability of the system being evolved by a unitary evolution.In other words, increasing the dissipation strength µ lowers the possibility of exponential growth which is evident in Fig. 6.

Pole structures of auto-correlation and spectral function
The two sets of Lanczos coefficients generically modify the behavior of the autocorrelation function and correspondingly its Fourier transform, known as the spectral function.It is interesting to investigate the pole structure of the auto-correlation function, similar to its closed system counterpart [4].In particular, our above analysis suggests that we might devise a generic form of an auto-correlation function.Let us assume a form of the auto-correlation function ) The pole of the auto-correlation function in the complex t plane.The black crosses indicate the poles without dissipation i.e., µ = 0 or the leading order term in (6.7), with the blue-shaded region as the region of analyticity.The red crosses indicate the poles with auto-correlation of the form of (6.6).In the weak coupling regime, the poles move away from each other along the y axis, thus resulting in an effect only a linear shift with magnitude µ/α 2 .
with α being some constant which is independent of µ.This form of the autocorrelation function looks similar to (4.5)-(4.6). 11We can easily see C(0, t) = sech(αt) reduces to the known closed system counterpart [4,15] and C(µ, 0) = 1.In other words, we can take this as a two-variable function and forget about its origin.Thus it is valid for any µ, not necessarily small.In fact, this function gives an asymptotic linear behavior of both a n = iαµ(2n + 1) ∼ iχµn and b n = αn of the form (6.1), without making any approximation of µ.Note that for µ = 0, the auto-correlation (6.5) reduces to C(0, t) = sech(αt), which is obtained for a n = 0 and b n = αn [4,15].To investigate the pole structure of (6.5), we set it to zero and find that the closest pole is located as The pole is not exactly at the imaginary t axis, rather is it shifted (see Fig. 7).
Then a reasonable question is to ask which kind of system has such a form of autocorrelation?Our answer is the dissipative SYK in the weak coupling regime, modeled by any generic random p-body Lindblad operators.Then our α dictates the system strength while µ encodes the dissipative strength.We have already found such results both analytically and numerically in previous sections.Only then, does the pole structure of (6.6) give the information of the operator growth.Thus, to connect with such growth, we expand (6.6) in the small µ regime as The o(µ 2 ) terms do contain both real and imaginary parts but they are not relevant to our discussion.In fact, (6.6) does suggest that the pole is not only squeezed by a factor of √ α 2 + µ 2 along the imaginary t axis but also shifted by a length of sinh −1 (µ/α)/( √ α 2 + µ 2 on the direction of negative x axis.The combined effect has a diagonal shift (see Fig. 7), squeezing the poles into the domain of the analyticity of (6.5) for µ = 0.However, in the weak dissipation regime, the squeezing is no longer required in the leading order and the closest pole structure still gives the growth of b n and thus the Krylov exponent.The effect of dissipation merely affects a linear shift of the pole.Of course, at zero dissipation, the poles are exactly located at t = ±iπ/(2α) [4].
The spectral function is given by the Fourier transform of the auto-correlation function, i.e., Taking the auto-correlation of the form (6.5), we find . (6.9)This is also a completely generic result, for any α, µ, and ω.However, we want to connect it to the operator growth and thus in the weak dissipation regime, Eq. (6.9) becomes The leading term in (6.10) is the known result for µ = 0 as for large ω , (6.12) which shows a long exponential tail in the large frequency regime [4].The subleading term in (6.11) depends on the ratio µ/α, and as long as µ ≪ α, the leading term dominates.Note that the leading term decays as exp(−#|ω|), while the subleading term decays as ω exp(−#|ω|).The overall decay still follows the leading behavior for large ω.Our analysis assumes a smooth behavior of the Lanczos coefficients unlike [16,29], where the different behavior of even and odd coefficients may lead to incorrect consequences on the spectral function. 12Thus, in the weak dissipation regime, the linear shift of the pole in the auto-correlation function is sublinear in µ, which reflects an ω exp(−#|ω|) decay in the sublinear term of the spectral function.This posits an alternate form of the operator growth hypothesis in open quantum systems.
7 Lindbladian dynamics: OTOC and q-complexity In this section, we derive an analytic expression for the out-of-time-order correlator (OTOC) for a time-evolved p-body fermionic operator and some general aspects of q-complexity.

Lindbladian dynamics for OTOC
We start with the OTOC of a single-body fermionic operator.The expression to be evaluated is the following where is the normalisation factor and ρ ∞ = 1 2 N /2 I.The overall factor of 1/2 arises because we use unnormalized ψ j , and the 1/N factor is introduced to account for averaging over the full set of Majorana fermions [2].The key idea is simple: the knowledge of the Krylov basis of ψ 1 allows us to write the time-evolved operator ψ 1 (t) in Krylov basis as follows In general, this would not give us a lot of analytical prowess.In the limit of large q and large N , however, things become tractable due to operator concentration (4.2).
Recall that the Krylov vectors are naturally orthonormal to each other.Using this we can evaluate the expression for the commutator {ψ j , ψ 1 (t)} as A general property of Krylov vectors is that basis vectors generated from a Hermitian Hamiltonian and initial operator are alternatingly Hermitian and anti-Hermitian.This allows us to evaluate the Hermitian conjugate of the above expression quite simply by noting that the combination i k O k is Hermitian.In other words, The trace operation under the sum is written as where we have used the fact that ρ ∞ = 1 2 N /2 I.The final task is to evaluate the trace operation in the RHS of the expression above, i.e., Note that the first term only contributes if k = l and (assuming that the Krylov vectors are properly normalized) contributes (−1) k 2 N /2 .The second term is the one that has to be evaluated carefully.We note that the trace of the product of fermion strings of different lengths vanishes.This implies that we have a non-zero contribution coming from k = l only.Therefore the full expression becomes We now look at the term ψ j O k in the second trace.It is straightforward to see that the following result is true Therefore the second trace term becomes where in the last step we have used the fact that c * i 1 i 2 ...is = (−1) kq/2 c i 1 i 2 ...is in order to ensure that i k O k is Hermitian.Combining the two terms we get the following expression Utilizing this, the trace operation under the full sum becomes We also note that The last piece of the ingredient is noting that since the Krylov basis vectors are normalized, it follows that . From these pieces of information, we can write the following expression for the OTOC ) .

(7.11)
There is a natural bound on this quantity, which follows from the fact that the sum of the coefficients in the expression satisfy In order evaluate the expression (7.11) further, we evaluate the following sum The first sum is simple.It involves performing the sum over the indices {i 1 , . . ., i s } over N − 1 values (i.e.excluding j).This evaluates to The second sum is non-trivial.To evaluate this we consider the cases where j ∈ {i 1 , . . ., i s }.This is broken up into two pieces, corresponding to j < s and j ≥ s.For j ≥ s, any i m can be chosen from i 1 , . . ., i s and set equal to j.Once the i m is chosen, the sum breaks into two pieces.
The remaining sum is Therefore in both these cases, the sums evaluate to the same value.Thus the full sum is Now, we note that under disorder averaging, all the coefficients |c i 1 i 2 ...is | 2 must be equal.Given that their sum is 2 s , each individual term is equal to Inserting (7.17) and (7.18) in (7.11), we obtain the following expression where in the final step we have inserted the expression for s in terms of k, noting that it is odd for even q.Therefore an analytic expression for the OTOC is given by We note that for the Lanczos coefficients endowed with the fairly generic form b 2 n = (1 − u 2 )n(n − 1 + η), a n = iu(2n + η), the Krylov basis wavefunctions φ n are given by (6.2).Using this in (7.20), we obtain the following expression (we kept some separation from u = 0 line to have the visual distinguishability) of a closed system.We choose η = 1 and q = 300.The behavior is similar to (6).
Fig. 9 shows the behavior of OTOC with different dissipative strengths.At late-times (t → ∞), the OTOC saturates to A useful estimate of the saturation timescale is the time at which the saturation value of the OTOC for a given u is equal to the u → 0 limit of the OTOC.It is straightforward to see that this timescale is given by This timescale is logarithmic in terms of the inverse dissipation strength.The analysis so far has been considered for the time-evolved initial operator ψ 1 .It is straightforward to see that exactly the same result holds for any single-body initial operator ψ i .However, the same results will not hold exactly for the general p-body initial operator of the form ψ i 1 ψ i 2 . . .ψ ip with i 1 < i 2 < ⋅ ⋅ ⋅ < i p .As we have discussed, operator concentration holds for any initial operator string.Therefore, most of the discussion will follow through with minor modifications (like replacing 1/ √ 2 in (7.2) by 1/2 p/2 ).One subtlety to be pointed out is that this c * i 1 i 2 ...is = (−1) (kq−p(p−1))/2 c i 1 i 2 ...is for a p-body initial operator.The normalisation factor is now This cancels the factor of 1/2 p/2 mentioned before.
The remainder of the calculation described above goes through in an identical fashion.The difference occurs at the final step of (7.19).The OTOC is there given by When s is odd (i.e. if p is odd), the term in the bracket becomes 2s/N , and for even s (i.e. for even p) the term becomes (2 − 2s)/N .However, note that the term 1/N just contributes a constant 1/N .Thus we are left with This concludes the discussion on OTOCs.

Aspects of q-complexity
We briefly discuss the formal results expected for q-complexity.It is a generalization of the operator growth probes like Krylov complexity, OTOC, and operator size.The notion is defined in detail in [4].We briefly describe the same below : • For a positive semi-definite superoperator Q, one can write it in its' own eigenbasis as follows • There exists a number K such that This condition ensures that the expectation value of Q does not change massively under one application of the Liouvillian.
• Similarly, there exists some number K ′ for an initial operator O such that This condition ensures that the initial operator has a low value of the complexity.
The q-complexity of an initial operator O is then defined as For our purposes, in the large −N large −q SYK model, it suffices to realize that any general superoperator in the space of Majorana fermions can be written as ..,is,j 1 ,...,jt q i 1 ,...,is j 1 ,...,jt |ψ i 1 . . .ψ is )(ψ j 1 . . .ψ jt | .(7.30) Due to operator concentration, one can write the time-evolved operator O(t) as follows where + p for a p-body initial operator.Using this expression, we find the time evolved q-complexity as follows The inner products in the above expression are evaluated by assuming that all the operators are properly normalized (or factors absorbed in ) and using the fact that strings of different lengths are orthogonal.This simplifies the expression to the following For q even, the string length d(k) can be either odd or even for a fixed p for all k.This is as far as one can go in general.However, using the disordered nature of our system, we can make the reasonable assumption that (with small variance and zero mean) under disorder averaging one can write ⟨c * i 1 ,...,i d(l . This greatly simplifies the expression for the (disorderaveraged) Q(t) further, leading to Further computation would require exact knowledge of the coefficients q i 1 ,...,i d(l) i 1 ,...,i d(l) .Specific choices of these coefficients lead us to the expressions for the K-complexity, OTOC, operator size, etc.For these three probes the exact coefficients are discussed in [4].

Summary and conclusions
In this paper, we performed a detailed study of operator growth through Krylov complexity in Lindbladian SYK, where the dissipation is modeled by various jump operators in the Markovian regime.In particular, we choose our system to be SYK and model the dissipation by p-body Lindblad operators.We analytically find the Lanczos coefficients which are numerically verified by implementing the bi-Lanczos algorithm, a suitable generalization of the Lanczos algorithm in open systems.We obtained a universal result of Krylov complexity, which initially grows exponentially and saturates at late times.Both the saturation time scale and saturation value appear to be generic and independent of the choice of the Lindblad operators, similar to what we obtained from OTOC.We also provide a plausible explanation of our results from the quantum measurement perspective.Based on these findings and analyzing the generic pole structure of auto-correlation and high-frequency behavior of the spectral function, we propose an operator growth hypothesis for generic open quantum systems, which suggests a broader notion of "dissipative quantum chaos".
Our approach opens a door to understanding a generic study of operator growth and chaos in non-Hermitian systems.Especially, since any Hermitian matrix can always be tri-diagonalizable and put into the form (3.8) with a n = 0, we wonder what could be a Hamiltonian structure of (3.8).In particular, for the Lindblad evolution, the density matrix evolves as [70] where H eff is known as the effective Hamiltonian which is non-Hermitian. 13This non-Hermitian Hamiltonian is constructed by the jump operators. 14It is important to note that the first term H is Hermitian while the second term is anti-Hermitian, which makes the overall Hamiltonian H eff non-Hermitian.Hence, we wonder whether any non-Hermitian Hamiltonian (including PT-symmetric [71,72]) can be cast into the structure of (3.8), which can be obtained by an efficient implementation of the bi-Lanczos algorithm.A limitation of our analysis is that we are completely blind to the physics with stronger coupling.A more broad analysis is to consider a generic coupling, not necessarily a small one.In fact, some preliminary analysis shows that the Lanczos coefficients become unstable at stronger coupling, which results from the fact that the Markovian approximation breaks down.It will be interesting to explore the non-Markovian regime where a generic coupling can be chosen.In fact, detailed knowledge about the environment (which might be another SYK) might lead to the study of the system in a purely scrambling or dissipative phase [42], resulting in an environment-induced phase transition.Finally, our bigger aim is to develop a systematic and coherent picture to understand dissipative quantum chaos in holographic duality.In particular, since the dual picture of generic p-body dissipator is unknown, it is interesting to see if our study in open SYK leads to the appreciation of a clearer picture of de Sitter space in quantum gravity.We hope to address this question in future studies.
Note added: During the final stages of this work, Ref. [67] appeared which deals with the behavior of the high-frequency regime of the auto-correlation.Their choice of auto-correlation function, specific types of system, and dissipation is fundamentally different from ours.Hence, we do not make any comparison with our results with Ref. [67]. where A slight renaming and rearranging of the indices gives us the following expression Clearly, this term is not ∝ ψ 1 .However, each term in the parenthesis is also a random variable with the same mean and the variance as V ij .Thus, under the condition that V k ij V k pq , (with i ≠ j and/or p ≠ q) vanishes upon averaging (for small V ), the second and third terms in the above expression vanish and we are left with Moreover, we have the result (A.9) In the large N limit, this concludes (A.5).
Now we consider a general string of length s = n(q − 2) + 1 and state the following proposition (4.12).
Proposition 3. The action of the dissipative part of the Lindbladian to a string of length s = n(q − 2) + 1 results in the following expression where ζ ∼ o(1) number and V is given by the ensemble average (4.11).
Proof: Here we discuss the general case for an operator string ψ i 1 ψ i 2 ⋯ψ is , where s = n(q − 2) + 1 is an odd number.For this, we note from the expression (A.1) that the only terms that will survive after averaging over multiple realizations are the ones with i = p and j = q.We only evaluate terms of this kind below.There are three such distinct cases.
Case 1, {i m , i n , i m , i n }: Here all four indices {i m , i n , i m , i n } ∈ {i 1 , i 2 , ⋯, i s } and the two pairs are identical.
Here we find For compactness, we denote the sum ∑ 1≤α 1 <α 2 <⋅⋅⋅<αp ans ∑ {α,p} in what follows.The dissipative part of the Lindbladian will be the central focus of our analysis.This is written as where we have used the following fact We emphasize that the coefficients V k α 1 α 2 ...αp are taken from a random complex Gaussian distribution with zero mean and variance ⟨|V k α 1 α 2 ...αp | 2 ⟩ = p!V 2 N p .This implies that for small V and large N , we can ignore terms of the kind ..βp with some indices different between the sets {α p } and {β p } since these will vanish upon averaging over a large number of realizations.We focus our attention to terms with α i = β i ∀ i.With these terms only, the average over Lindbladian may be written as The terms in the squared parenthesis need to be treated with care.We split these into two terms where we have used the expression for O which follows from operator concentration.Firstly, we note that The evaluation of G 1 is more involved.For this, we note that any α k can be written as Systematically moving each α i across ψ i 1 ψ i 2 ⋯ψ is from the right to left (or vice-versa) we pick up a phase of (−1) p−i by the time it reaches ψ i 1 ψ i 2 ⋯ψ is .At which point it crosses over to the other side picking the phase (−1) s+ζ s i .Then it combines with the corresponding operator on the other side of ψ i 1 ψ i 2 ⋯ψ is to give a factor of 1  2 I. Repeating this process for all the p fermions, we find the following result Combining G 1 and G 2 , we get the following With this, one can write the averaged dissipative Lindbladian as follows One feature that we note from the onset is that if none of the α i indices lie in {i 1 , i 2 , . . ., i 2 }, then we get zero contribution since all the ζ s l are vanishing.This was also observed explicitly in the 2-body jump operator case studied in the previous section.Additionally, it is also evident that there will be a non-zero contribution only when an odd number of the indices {α i } lie in {i 1 , . . ., i s }.This was also observed in the 2-body jump operator case.
The objective now is to perform the combinatorial sum It is useful to consider this sum in pieces.This sum tells us to put 2 for every case where 1 ≤ α 1 < ⋯ < α p ≤ N and then subtract 2 for each case where ∀α i ∉ {i 1 , . . ., i s }.
Then the cases where an even number of α i ∈ {i 1 , . . ., i s } also have to correspond to a subtraction of 2. We denote the L.H.S of the above summation as 1 2 (S 1 − S 2 − S 3 ).To evaluate the first step, we note the following sum The next step is to subtract out the cases where ∀α i ∉ {i 1 , . . ., i s }.To do this, it is enough to note that the sum will be the same as (B.16), just with N replaced by N − s.Therefore this term is . (B.17) Let us now look at the term S 1 − S 2 .In the large N limit, this terms is given by The final step is to evaluate the sum where an even number of α i are taken from {i 1 , . . ., i s }.Since the numbering of the fermions is really arbitrary, it suffices to compute the expression for the arrangement {i 1 , . . ., i s } = {1, 2, . . ., s}.Any other arrangement of the indices will give the same result.Since we evaluate the sum for the case where {i 1 , . . ., i s } are sequential, the combinatorial factors will arise from the choice of an even number of {i 1 , . . ., i s }.There will be no combinatorial factors from assigning these indices to {α 1 , . . ., α p }, since these will necessarily be assigned in an increasing order starting from α 1 .
The number of ways of selecting 2k number of indices from the s indices available to us is ( s 2k ).Once these 2k elements are chosen and assigned to α 1 , . . ., α 2k in ascending order, the sum can be performed over the remaining p − 2k summation indices spanning over N −s values.Note that 2k is limited by the s or s−1 (depending on s odd or even) or p (p − 1 if p is odd).This simply amounts to replacing N by N − s and p by p − 2k in (B.16).The resulting summation is ))ψ i 1 ψ i 2 ⋯ψ is . (B.20) In the large N approximation, we ignore S 3 as it is easy to see that the leading order contribution of S k is O(N p−2k ), which is at least one order less than N p−1 which is the leading order contribution from S 1 − S 2 .Hence, It is easy to see that this reduces to the leading order contribution (A.17) for p = 2. Hence, the "operator size concentration" (4.2) leads to (4.14).

C General Initial operator
It has been shown that operator concentration is a property of the large−N , large q SYK model [22].Using half-melon diagrams to represent the operators generated by subsequent application of the closed SYK Lindbladian, starting with initial operator ψ 1 , it was shown that the Lanczos coefficients are given by b 1 = J √ 2/q, b n>1 = J √ n(n − 1).It is straightforward to derive the same results for a general p−body initial operator ψ i 1 ψ i 2 . . .ψ ip .Before presenting the arguments, let us briefly review the salient features of the construction in [22].
The SYK Louivillian L H can be split into two parts, each corresponding to a forward and backward movement in the Krylov basis respectively Considering an operator basis generated by the action of L + on the initial operator, it is possible to demonstrate that the basis is orthonormal.This is due to the fact that L + generates Majorana strings of length larger than the one on which it acts.The largest such string is one that arises when a Majorana fermion (in the operator on which L + is acting) is replaced by a (q − 1) body string of fermions.A diagrammatic approach is discussed in detail in [22].Subsequent levels of L n + ψ 1 give rise to a rapidly increasing number of diagrams, each of the same length.The red curve is the initial operator and the further black curves represent the subsequent additional operators generated.The coefficients of each diagram are essentially the number of ways it can be constructed out of its parent diagram in the previous level.Careful counting results in the following identity Using this, one can argue that the Krylov basis up to normalization is given by And the usual action of the Liouvillian on the basis element follows accordingly (for n > 1) with the edge cases n = 0, 1 represented as The edge cases have to be calculated by explicit calculation as follows.We first consider the case n = 0, remember we start with the normalized initial operator The commutator in the sum [ψ i 1 . . .ψ iq , ψ 1 . . .ψ p ] evaluates to (1−(−1) l )ψ i 1 . . .ψ iq ψ 1 . . .ψ p where l is the number of the indices in i 1 , i 2 , . . ., i q that coincide in some index from 1, 2, . . ., p.This implies there is a non-zero contribution only when l is odd.The commutator is then denoted as Using the definition of the variance ⟨|J i 1 i 2 ...iq | 2 ⟩ = (q−1)!J 2 N q−1 , and the redefinition 2 1−q J 2 = J 2 q corresponding to the large−q limit we can simplify the summation in b 2 1 under disorder averaging to the following expression b 2 1 = 2 2−q (q − 1)!2 q−1 N q−1 J 2 q ∑ 1≤{i,l}≤N In the large N limit, we obtain the following expression b 2 1 = 2 2−q (q − 1)!2 q−1 N q−1 J 2 q ∑ 1≤{i,l}≤N 1 → N →∞ 2 Γ(q) N q−1 qN q−1 p Γ(q + 1) (C.28) (a) Compute: |r j ⟫ = L † o |p j ⟫ and |s j ⟫ = L o |q j ⟫.(b) Redefine the vectors: |r j ⟫ ∶= |r j ⟫ − b j−1 |p j−1 ⟫ and |s j ⟫ ∶= |s j ⟫ − c * j−1 |q j−1 ⟫.(c) Evaluate the inner product: a j = ⟪q j |r j ⟫.(d) Again, redefine the vectors: |r j ⟫ ∶= |r j ⟫ − a j |p j ⟫ and |s j ⟫ ∶= |s j ⟫ − a * j |q j ⟫.(e) Evaluate the inner product: ω j = ⟪r j |s j ⟫.(f) Evaluate the norm: c j = √ |ω j | and b j = ω * j /c j .(g) If c j+1 ≠ 0, then define the vectors:

Figure 1 .
Figure 1.The operator growth in dissipative systems (left) can be mapped to a model of the non-Hermitian Krylov chain (right).The hopping amplitudes from n-th side to (n + 1)-th and (n − 1)-th sites are b n+1 and c n respectively, while a n+1 gives the amplitude of staying at site n.Here O 0 indicates the initial operator and the red arrows indicate (increasing length indicates stronger dissipation) the dissipation which affects all sites.
Behavior of |a n |.
Behavior of b n .

Figure 2 .
Figure 2. Behavior of the (a) diagonal coefficients |a n | and the off-diagonal coefficients b n with different dissipative strength for the SYK 4 model, with linear dissipators.The dotted line in (a) is given by (5.1).Our initial operator is O 0 = √ 2ψ 1 , and the number of fermions, N = 20.Here we have taken 50 Hamiltonian realizations.
Behavior of |a n |.
Behavior of b n .

Figure 3 .
Figure 3. Behavior of the (a) diagonal coefficients |a n | and (b) the off-diagonal coefficients b n with different system sizes for SYK 4 , with linear dissipators.The insets show the linear dependence (fitted) of the saturation values of |a n | and b n .Our initial operator is O 0 = √ 2ψ 1 , and the dissipative strength is fixed at λ = 0.01.Here we have taken 50 Hamiltonian realizations.
and Fig. 5b, we separately show the behavior of the slope m(|a n |) and the saturation value |a (sat) n | with the number of Lindblad operators.This increase is linear in either case, i.e., m(|a n |) ∝ M , and |a Behavior of |a n |.
Behavior of b n .

Figure 4 .
Figure 4. Behavior of the (a) diagonal coefficients |a n | and (b) the off-diagonal coefficients b n with the different number of jump operators M for SYK 4 .Notice that the b n s exactly overlap for all values of M .Our initial operator is O 0 = √ 2ψ 1 , system size is N = 16 and the dissipative strength is fixed at V = 0.02.The random Lindblad operators are taken over 30 averages.

Figure 5 .
Figure 5. Behavior of the (a) slope of the diagonal co-efficient |a n | and (b) the saturation of the diagonal coefficients with the different number of jump operators.Our initial operator is O 0 = √ 2ψ 1 , system size is N = 16 and the dissipative strength is fixed at V = 0.02.The values are obtained after averaging over 30 disordered realizations.(c) Dependence of c V with V .The numerical fitting gives κ = 0.0780 and β = 2.0046 according to (5.6).

Figure 6 .
Figure 6.The behavior of Krylov complexity (6.3) for different dissipation strengths.The gray dashed line indicates the behavior of ∼ e 2t (we kept some separation from u = 0 line to have the visual distinguishability) of a closed system.We choose η = 1.

Figure 9 .
Figure 9.The behavior of OTOC(7.21)  for different dissipation strengths.The gray dashed line indicates the behavior of ∼ e 2t (we kept some separation from u = 0 line to have the visual distinguishability) of a closed system.We choose η = 1 and q = 300.The behavior is similar to(6).