Universal equilibration dynamics of the Sachdev-Ye-Kitaev model

Equilibrium quantum many-body systems in the vicinity of phase transitions generically manifest universality. In contrast, limited knowledge has been gained on possible universal characteristics in the non-equilibrium evolution of systems in quantum critical phases. In this context, universality is generically attributed to the insensitivity of observables to the microscopic system parameters and initial conditions. Here, we present such a universal feature in the equilibration dynamics of the Sachdev-Ye-Kitaev (SYK) Hamiltonian -- a paradigmatic system of disordered, all-to-all interacting fermions that has been designed as a phenomenological description of quantum critical regions. We drive the system far away from equilibrium by performing a global quench, and track how its ensemble average relaxes to a steady state. Employing state-of-the-art numerical simulations for the exact evolution, we reveal that the disorder-averaged evolution of few-body observables, including the quantum Fisher information and low-order moments of local operators, exhibit within numerical resolution a universal equilibration process. Under a straightforward rescaling, data that correspond to different initial states collapse onto a universal curve, which can be well approximated by a Gaussian throughout large parts of the evolution. To reveal the physics behind this process, we formulate a general theoretical framework based on the Novikov--Furutsu theorem. This framework extracts the disorder-averaged dynamics of a many-body system as an effective dissipative evolution, and can have applications beyond this work. The exact non-Markovian evolution of the SYK ensemble is very well captured by Bourret--Markov approximations, which contrary to common lore become justified thanks to the extreme chaoticity of the system, and universality is revealed in a spectral analysis of the corresponding Liouvillian.


Introduction
Whether and how a perturbed system equilibrates have been fundamental issues of statistical mechanics since the laying of its foundation. One question that has intrigued researchers for almost a century is the thermalization of an isolated quantum system under its unitary evolution [1][2][3][4]. In the last two decades, this process has experienced a revitalized surge of interest, thanks to experimental breakthroughs in realizations of synthetic many-body systems [5][6][7][8][9][10][11]. An unprecedented control over system parameters now enables laboratory investigations using quantum systems in almost ideal, isolated conditions [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. On the theory side, a main obstacle for arriving at a unified understanding of out-of-equilibrium quantum dynamics is the absence of a universal principle that would be as general as the minimization of free energy for equilibrium phase transitions [27,28]. Nevertheless, powerful frameworks have been developed to explain the thermalization of a quantum system, perhaps the most successful being the eigenstate thermalization hypothesis (ETH) [29][30][31][32]. According to the ETH, for quantum chaotic systems, particularly for ergodic systems, thermalization and equilibration are tantamount, since under quantum chaotic dynamics a perturbed system equilibrates to a state that for local observables is indistinguishable from a Gibbs thermal state. A well accepted mechanism for thermalization in an isolated quantum system is that initially localized information is distributed among the system's degrees of freedom, and thus becomes irretrievable through any local operation at later times. This process, referred to as scrambling [33,34], occurs in quantum many-body lattice models [35][36][37], conformal field theories [38], and black holes alike [39,40], and has been experimentally probed in different physical systems [41][42][43].
In this paper, we identify a universal equilibration in quench dynamics of the complex SYK model, revealed in state-of-the-art numerical calculations for the exact dynamics and reproduced analytically through a master equation. In particular, employing a highly optimized exact diagonalization method for systems comprising up to N = 20 complex fermionic modes, we study how a system initialized in an eigenstate of some other Hamiltonian equilibrates to a steady state following a sudden global quench into the SYK model. For a broad variety of few-body observables, including multipartite entanglement as given by the quantum Fisher information (QFI), the disorderaveraged evolution collapses onto a single curve after a simple amplitude rescaling, independent of (generic) initial states. Over vast stretches of the dynamical evolution, this universal curve is well approximated by a Gaussian, with a fast decay on the order of the time-scales of leading-order processes. Such a universality over the entire evolution goes significantly beyond what is observed in conventional equilibrating systems, where universal behavior independent of initial conditions can only be expected once the system reaches a fixed point of the dynamics (the final steady state [3,4,32] or a non-thermal fixed point [25,[88][89][90][91][92][93]). Thus, our findings may stimulate future investigations in non-equilibrium quantum many-body dynamics in order to identify similar universal dynamics in other models.
We substantiate our numerical findings by devising a Lindblad master equation (ME) that describes the Hamiltonian disorder average as an effective nonunitary time evolution. In this formalism, the unitary but disordered closed-system dynamics generated by the SYK model is mapped to one of a clean but dissipative system. A detailed prescription for Hamiltonian disorder averaging has been introduced by Kropf et al. based on a matrix formalism [94][95][96]. Here, we present an alternative, mathematically elegant route to Gaussian disorder averaging based on the Novikov-Furutsu theorem. In earlier works, this theorem has been applied in the context of averaging noise with finite correlation times in, for instance, quantum walks subjected to pure dephasing noise [97,98], stochastic Schrödinger equations [99], or proposals for simulating dissipation via noisy unitary dynamics [100][101][102][103]. We exploit that framework by formally promoting the quenched disorder to noise with infinite correlation time, permitting a fruitful application to a generic system with Gaussian disorder. To render the ensuing exact equations tractable, we employ decorrelation and Markovian approximations. In contrast to standard lore, which seemingly preempts their use for the infinite correlation time of quenched disorder [104][105][106][107][108][109][110], these approximations capture the true quantum dynamics well, thanks to the extreme chaoticity of the SYK model. The resulting master equation successfully describes the super-exponential aspect of the equilibration process as well as the equilibrated steady state. Even more, a spectral analysis of the associated Liouvillian provides valuable insights into the universal dynamics of few-body observables. As these results show, the master equation for the ensemble averaged state elegantly describes complex features of disordered quantum dynamics, and may thus constitute a powerful tool beyond the immediate context of this work.
The rest of this manuscript is organized as follows. In Sec. 2, we elaborate on a global quench protocol, employed to induce the equilibration dynamics. Then, in Sec. 3, we discuss the universality and superexponential decay observed in the disorder ensemble averaged dynamics of the QFI. Sec. 4 presents the formalism used to obtain the master equation, and further illustrates that the latter captures the salient features of the equilibration process. These sections constitute the main results of our study. To provide a further in-depth analysis, we extend our study in Sec. 5 to the dynamics of operator moments, and explain the observed universality in the equilibration process based on a spectral analysis of the Liouvillian. In Sec. 6, we conclude the key findings of our study, and emphasize possible extensions and potential applications of the presented formalism. The main article is complemented by appendices that provide further details on the model, master equation formalism, as well as additional numerical results.

Quench protocol
We are interested in the disorder-averaged out-ofequilibrium dynamics generated by the SYK q model for spinless complex fermions, following the quench protocol sketched in Fig. 1a. The family of SYK q models is characterized by (q/2)-body random all-toall interactions, where q/2 ∈ Z + . In this paper, we present in detail the dynamics of the SYK 4 model, but we stress that the salient features of this equilibration process are qualitatively generic to other values of q [111]. The Hamiltonian of an instance of the ensemble reads (see App. A for details) [48,112] There are N spinless fermionic modes with creation, annihilation, and occupation number operators that satisfy canonical anticommutation relations and are denoted, respectively, byĉ † i ,ĉ i , andn i . The interaction strengths {J i1i2;j1j2 } are complex Gaussian random variables, and we denote the disorder average over realizations by E[...].
As initial states |ψ(0) , prepared at time t < 0, we consider for convenience ground states of the one-dimensional spinful Fermi-Hubbard (FH) model [113,114] given by the Hamilto- =1n ,↑n ,↓ . In the FH model, physical modes are given by N/2 spatial lattice sites, = 1, . . . , N/2, and additional spin degrees of freedom (represented by the arrows) not present in the SYK model. As the family of SYK q models consists of zero-dimensional models (due to site-independent all-to-all interactions), the mapping of FH modes to the SYK modes is arbitrary. We choose here { , ↑} ↔ i = 2 and { , ↓} ↔ i = 2 − 1. The FH ground states are chosen as a representative case of initial states whose operator expectation values differ significantly from the steady-state values.
In our numerics, we employ various choices of the boundary condition, 1 ratio between onsite interaction and hopping strengths U/J , total fermion number 1 ≤ Q ≤ N , and total magnetization (we focus mostly on the case of half-filling, Q = N/2, and zero magnetization, where there are N/4 fermions in each spin sector). We emphasize that the above choice of the initial Hamiltonian is only for convenience of preparing initial states that cover a range of parameters in a strongly-correlated system. Choosing other generic initial states does not modify our findings.
Once the system is prepared in the initial state |ψ(0) , we perform a global quench at t = 0 to the SYK 4 model and track the state's subsequent unitary time evolution (here and throughout we set = 1) We average the time series over multiple disorder realizations in order to filter out the salient, realizationindependent features of the equilibration dynamics.
To test the generality of our findings, we study a variety of observables, in particular the QFI (see next section) and higher-order correlators (see Sec. 5 and App. B) corresponding to the staggered magnetization, as well as other few-body (4-local) operators and a non-diagonal generator for the QFI (see App. C). As a main result of our study, the disorder-averaged evolution of the considered few-body observables 1 We performed computations for periodic, anti-periodic, and open boundary conditions. In this article, we present the dynamics of initial states obtained with periodic and antiperiodic boundary conditions when N/4 is odd and even, respectively.
shows universality under the following rescaling where f (t) represents the long-time average of the function f (t) computed over a time window starting at t 0 and with duration T , i.e., Unless explicitly stated, we consider Jt 0 = 50 and J(t 0 + T ) = 100 for the exact diagonalization results presented in this paper. It is to be noted that the SYK model is only parameterized by the variance of the random all-to-all interaction strengths [see Eq. (27)]. Therefore, the universality in the evolution of observables can be probed as the insensitivity to the initial conditions, which becomes evident under the rescaling in Eq. (3), as is discussed throughout this manuscript.

Universal super-exponential equilibration dynamics
To illustrate the universal equilibration dynamics, we present here results for the QFI evolved under the SYK 4 Hamiltonian. The QFI is an observable of central relevance in quantum sensing [115][116][117], which can witness multipartite entanglement in quantum many-body systems at zero and finite temperatures [118][119][120][121][122]. Interestingly, like the out-of-timeorder correlators [123,124], this quantum information theoretic measure can distinguish a pure eigenstate of an ETH-obeying Hamiltonian from the corresponding Gibbs thermal state [125].
In the present context of pure states, the QFI with respect to an observableÔ is simply proportional to its variance (5) In this section, we consider the staggered magnetization, which in the FH model is defined asÔ SM = N/2 =1 (−1) (n ↓ −n ↑ ), and which in the SYK model translates tô Theκ i denote 2-local operators which we use to construct the 4-local operators discussed in App. C. The time evolution of the disorder-averaged QFI, E [F Q ], computed with respect to the operatorR is shown in Fig. 1b. The considered initial states are the symmetry unbroken ground states of the FH model for U/J = 0, 2, 4, 6, 8, and 10, respectively (dark to light shading). These values include a non-interacting initial system for U/J = 0, as well as strongly interacting systems at larger values of U/J . As a result, the initial states are characterized by a varying amount of multipartite entanglement that is witnessed by the QFI. 2 At short times, the system still retains memory of the initial state. However, the completely disordered all-to-all interactions of the SYK model lead to a quick loss of this memory, and already at a time of about Jt ≈ 4 the QFI equilibrates to a steady state value that is independent of the initial state. This rapid equilibration is reminiscent of the fast scrambling characteristic of the model, and also bears similarities to the relaxation curves derived in Ref. [126] for the out-of-equilibrium dynamics of isolated quantum many-body systems. There it is shown that rapid, non-exponential equilibration is expected for, in a random matrix sense, typical Hamiltonians and observables.
The attained steady state value matches with the one obtained from the infinite temperature Gibbs state (horizontal dashed black line in Fig. 1b) where 1 is the identity operator, Z is the partition function and D is the Hilbert space dimension determined by N and Q. This finding suggests that the overlaps between a generic initial state and energy eigenstates of the SYK 4 Hamiltonian are uniformly distributed over the spectrum. This is substantiated by computing the Kullback-Leibler divergence, D KL (P (E) Q(E)), between the uniform distribution Q(E) = 1/D and the initial states' amplitude distribution P (E) = | ψ(0)|E | 2 with respect to the energy basis {|E }. 3 We note that even though the steady stateρ ∞ has vanishing QFI (see for instance Ref. [115]), it is nevertheless an interesting question of how the system reaches that point, starting from initial states with different amounts of quantum correlations. Indeed, we find that this equilibration dynamics is universal within numerical precision, as one can expose by rescaling E [F Q ] according to Eq. (3). The dynamics of the rescaled curves are shown in Fig. 1c. All the curves collapse throughout the dynamics, independent of the initial state. In addition, the universal curve can be well approximated by a Gaussian (red dashed curve in Fig. 1c), with a fast decay constant of 2 We note that the scaling of the QFI with system size N depends on the amount of multipartite entanglement of the state [117]: For separable states states, F Q ≤ N , whilst for genuinely N -partite entangled states F Q ≤ N 2 .  We note that within the context of random matrix theory it is well known that a Gaussian temporal evolution can occur in the survival probability. For instance, this happens for quantum quench dynamics under Wigner random banded matrices and two-body random ensembles (see for instance Ref. [127], and Review [128] and references therein), where the latter is a specific case of embedded random matrix ensembles [129]. The same behavior has also been studied for a generic disordered interacting spin model [130]. In these cases, the survival probability initially decays as a Gaussian, followed by a regime in which it shows oscillations with a power-law envelope [130]. The oscillatory behavior is also seen in the evolution of the spectral form factor of SYK models [57,63]. Our numerics illustrate the above features for the survival probability (see Fig. 11). Indeed, for q = 4 the SYK model can be interpreted as a two-body random ensemble (albeit without a mean-field contribution), and more generally as an embedded random matrix ensemble for q ≥ 2 [see Eq. (26)], and as such the above features are expected. In contrast to this many-body observable, for the few-body observables considered in this study, we find the Gaussian decay followed by a marginal domain in which an indication of a powerlaw tail is obtained, but which appears to diminish with system size (see Fig. 5).

Dissipative ensemble dynamics
In this section, we present the key results of the opensystem formalism, which we use to understand the main features observed in the average dynamics of the unitary ensemble discussed in the previous section. Due to its generality for treating disorder ensemble averages, this formalism is of interest also independent of the application to the present scenario. The derivation presented here is formally written for time-dependent random processes. We stress, however, that we allow for time-dependent stochastic processes in Eq. (8) only to make our formalism applicable to more general Hamiltonians. The resulting evolution equations apply also in the limit of static processes, such as those defining the SYK model. The approach presented here has the advantage of treating quenched disorder and temporally fluctuating noise on equal footing, enabling an application to a large variety of settings. The interested reader may find further details and the explicit derivation for static processes in App. D.
Our aim is to directly study the dynamics of the ensemble's density matrixρ(t) ≡ E [ρ(t)] via the ensemble averaged von Neumann equation (EAVNE), where the time-evolution of each stateρ(t) is generated by a realization of the general closed system Hamiltonian Here,Ĥ 0 (t) is a disorder-free contribution, which in general can be time dependent, whereas the termŝ capture the dynamics due to disorder or noise. The index α is used to distinguish different subsets of Hermitian operatorsĥ (α) lα , and the operators within a subset are labeled by the (multi-)index l α . 4 In particular, upon rewriting the SYK 4 Hamiltonian in the generic form of Eq. (9), we identify three operator subsets, as shown in Sec. 5.2. We assume the functions ξ lα,kα (t, t ) is constant with respect to time. To keep the formalism general, we will specialize to this case only further below.
Consider the EAVNE generated by averaging over multiple realizations of the Hamiltonian in Eq. (8), (10) To proceed, one needs to handle the correlations E ξ The simplicity of our framework rests upon use of the Novikov-Furutsu theorem [131][132][133][134], which provides an exact expression of these correlations in terms of F . (11) An explicit expression for the functional derivative can be obtained from the integrated von Neumann equation. Formally, this yields an infinite series in which the nth term (n ≥ 1) contains n − 1 time integrals over n nested commutators, from which the exact functional derivative can be obtained in principle. For a systematic study of the role of the higher order terms, we refer the reader to our follow-up work Ref. [135]. Here, we retain only the lowest order contribution (n = 1), which reads where the step-function Θ arises from causality. Substituting Eqs. (11) and (12) into Eq. (10), we obtain the evolution equation This evolution equation is not exact, due to our use of the approximate functional derivative given by Eq. (12). Such a leading order truncation amounts to the well-known decorrelation assumption 5 typically made in the analysis of stochastic evolution equations [136,137]. The remaining time integral in Eq. (13) is known as a Bourret integral [138,139]. While the decorrelation assumption becomes exact in the limit of white noise, for non-Markovian noise it corresponds to an expansion controlled by the noise correlation time [140]. One may then wonder what justifies this approximation (see also Fig. 2) for the present disordered system, which has an infinite correlation time. The reason can be attributed to the chaoticity of the SYK 4 model, rigorous proof of which remains an open question for future investigations. Each term of the Hamiltonian in Eq. (9) can be thought of as an independent noise process governing the evolution of the density operator. Then, in the presence of a large number of such processes-as in the SYK 4 modelone may expect the correlations between the density operator and any individual process to be strongly suppressed. Viewed differently, the decorrelation assumption can be seen as a linear-response approximation [104,105], i.e., the response of the stateρ at time t towards a perturbation with ξ lα at an earlier time t 1 is taken into account only to linear order. In the context of Kubo's celebrated linear response theory [106], it is well known that averaging effects due to chaos lead to a superb success much beyond the regime of applicability predicted by naïve estimates [104,[107][108][109][110]. In the present context, the observed success of the linear approximation can be seen as a manifestation of the strong effects of quantum chaos in the SYK 4 model. 5 This is readily seen within the interaction picture generated byĤ 0 (t), where the truncated functional derivative now acts on the transformed functional given by the interaction picture density matrix.
The master equation as given by Eq. (13) is still rather unwieldy, since it is not local in time. We thus perform a Markov approximation leading us to a Lindblad master equation in non-diagonal form [141] governed by a time-dependent Liouvillian superoperator (14) with Hermitian dissipator in which we have defined Equations (14)- (16) form the final evolution equations of this section. They are valid under a Bourret-Markov approximation for the generic Hamiltonians of Eq. (8) with disorder and/or noise contributions. We reiterate that these effective evolution equations do not require the presence of noise (see App. D, fourth paragraph), and that the master equation is a result of averaging over an ensemble of disorder realizations. Whilst each individual disorder realization evolves unitarily, the ensemble evolves like an open system, with a dynamics that is approximately generated by the Liouvillian L(t). The coherent and dissipative processes that constitute L(t) can be read-off immediately from the system's Hamiltonian. The corresponding dissipation rates are entirely determined by the disorder statistics F Similarly, while each realization preserves the purity of the initial state, the Hermitian jump operators of the master equation generate ensemble dynamics that are purity decreasing [142,143], and thus drive the ensemble from a pure to a mixed state. In particular, for the SYK 4 model with large enough N , the ensemble equilibrates to the infinite-temperature stateρ ∞ given in Eq. (7) [143,144]. Note, however, that whilst the Hermitian jump operators of Eq. (14) and (15) ensure thatρ ∞ is a steady state of the Liouvillian dynamics, for an arbitrary Hamiltonian as given by Eq. (8), the steady-state of L(t) need not be unique in general [145,146].
Besides clarifying the nature of the steady-state, the Liouvillian dynamics also reproduce the rapid, super-exponential equilibration of the QFI observed in the SYK 4 quench dynamics of Sec. 2: In the above equations, the case of static disorder is captured by a "noise" correlation that is constant in time, so that 2f kα . Thus, under the Bourret-Markov approximation the SYK 4 model (or indeed any SYK q model) is governed by dissipation rates that grow linearly in time. Consequently, one can factor the Liouvillian as L(t) = 2tD, which naturally yields the super-exponential time-evolutioñ where T denotes time-ordering. Figure 2a shows simulations of the QFI evolution generated by the above master equation. In general, the ME as developed here may not be used to study the ensemble average of the QFI, due to the term ψ(t)|Ô |ψ(t) 2 [see Eq. (5)] that is non-linear in the density matrix. However, for the initial zero magnetization states considered here, exact numerics show that the first moment fluctuates around zero, such that in this case one may use the ME to approximate the dynamics of the QFI. Indeed, the agreement with the ensemble averaged ED results is striking, and we emphasize that no fit parameter has been used (nor is one available in the above formalism) to achieve this agreement. The Liouvillian dynamics also reproduce the universality of the QFI with respect to different initial states, as shown in Fig. 2b. A discrepancy at intermediate times can be attributed to the decorrelation and Markov approximations of the master equation. These approximations can be expected to hold especially at early times, as can be seen from a shorttime expansion, and at late times when the system enters a steady state (a unique steady state is the one which is reached independent of the precise trajectory, so correlations to the exact noise process and memory about previous times can be expected to be negligible). Indeed, the ME provides excellent agreement at early and late times, and successfully captures the overall trend of the dynamics even at intermediate times.
Even more, the Liouvillian formalism allows us to study the origin of the universal dynamics. In Sec. 5.3, we analyze how different states and observables decompose over the eigenspaces of L(t), showing how these distributions conspire in the presented scenario to select a single time-scale, universal across different initial states.
As a final remark, the fact that disorder induces dephasing between members of an ensemble-and thus leads to effective open-system evolution equations even in the absence of a heat bath-is long known; in particular for special single-body cases such as classical disordered dipoles and harmonic oscillators [140] or single atoms coupled to a photon field [147]. Here, our aim was the derivation of a general framework for quantum many-body systems. Such a platform for the evolution of disorder-averaged density operators has been rigorously developed previously [96], based on a matrix formalism [94,95]. Our derivation based on the Novikov-Furutsu theorem is simpler but nevertheless general, as the assumptions we made are not fundamental, but rather of practical nature: The (a) For each initial state, the ME simulation reproduces the dynamics of the exact numerics very well. There is a discrepancy at intermediate times due to non-Markovian effects and higher-order correlations, not captured by the approximate ME (inset). (b) The ME reproduces the collapse to a universal curve under the rescaling defined in Eq. (3), without any free fit parameter. The slight spread in the rescaled ED curves is due to statistical fluctuations, which for the QFI are suppressed for larger system sizes, as can be seen from a comparison with Fig. 1c, which shows QFI dynamics computed via ED for N = 16.
Novikov-Furutsu formalism can be extended to non-Gaussian stochastic processes [134,148], the decorrelation assumption may be lifted in favor of an infinite series of terms in Eq. (12) [105,140], and the Markov approximation is-at least on a formal levelnot necessary in the derivation of the evolution equations. The present framework has the additional feature that disorder and noise processes can be treated on equal footing, within the same master equation, without any further complications of the formalism. That property enables us to tap into the vast literature employing the Novikov-Furutsu theorem in the context of noise with finite correlation time (see, e.g., Refs. [97,[100][101][102]).
To summarize this section, the Novikov-Furutsu theorem enables us to derive a master equation, for the ensemble average, that provides a general framework for disorder-averaged quantum many-body systems. For the case of the SYK 4 model, it yields a series of analytic insights into the out-of-equilibrium dynamics, such as the steady state reached at late times, the approximately Gaussian decay, and the universal behavior across different initial states.

Universal evolution of operator moments under SYK 4
To corroborate the generality of the above findings, in this section we consider the kth moment, M k (t) = ψ(t)|Ô k |ψ(t) , of the operatorÔ =R defined in Eq. (6) (in App. C, we report analogous results for 4-local operators and QFI computed with respect to a non-diagonal operatorT , defined in Eqs. (29) and (30), respectively). We start by presenting the corresponding numerical ED results which, as in the QFI case, display universal behavior in the disorderaveraged time series. Then, we use the Lindblad ME derived in Eq. (14) to further illuminate the universal dynamics and the salient features of the exact evolution. Finally, we present a spectral analysis of the corresponding Liouvillian, which explains the universality (initial state independence) within the Bourret-Markov approximation.

Numerical results from exact diagonalization
With respect to the symmetry unbroken FH ground states, the expectation values of all the odd moments of the staggered magnetization operatorR are zero. Their ensemble averaged expectation values continue to show negligibly small fluctuations around zero dur-ing time evolution under the SYK 4 Hamiltonian. 6 In contrast, the even moments exhibit the same superexponential universal equilibration behavior as the QFI. This is illustrated in Fig. 3 for k = 2, 4, and 6 (higher-order moments are presented in Fig. 9). For visualization purposes, in the left column the expectation values of the operator moments are normalized to values ≤ O(1) by an empirical factor The super-exponential approach to equilibrium is clearly visible in this data. As is evident from the right columns of Figs. 3 and 9, a rescaling using Eq. (3) collapses the even moments evolved from different initial states onto a single curve. During most of the evolution, this collapsed curve can be well approximated by a Gaussian. In the transient regime, curves corresponding to different initial states for even k ≥ 4 do show small deviations, an effect that becomes more prominent for larger N . In other words, while the universality found for k = 2 is very robust, for larger k it becomes approximate in an intermediate time window. Below, we describe this feature in detail via a spectral analysis of the Liouvillian. This finding also suggests that universality is more precise for few-body operators, as is further corroborated by comparison with the global many-body observable of the survival probability, see Fig. 11. However, we do not exclude the possibility of a suitably constructed, highly non-local observable exhibiting universality with respect to Eq. (3).
An interesting feature of the different moments is that their respective curves shift towards earlier times with increasing order k. That is, the higher the moment order is, the faster is the equilibration. Accordingly, Gaussian fits to the universal curves for k = 2, 4, 6 yield the decreasing decay times τ = 1.52, 1.42, 1.35, respectively, see Fig. 3. To illustrate this effect further, we show the rescaled curves for k = 2, 4, ..., 12 in Fig. 10, where it appears the curves converge with sufficiently high order. Again, this trend can be explained based on a spectral analysis of the Liouvillian as a function of moment order (see Sec. 5.3 and inset of Fig. 10 in App. B for details).
To study the finite-size dependency of the rescaled universal curves, we consider in Fig. 4 the representative case k = 2 of the operator moments. We consider the exact evolution for systems consisting of 6 One can prove formally that odd moments ofR vanish as follows. The considered initial states, having zero magnetization, have spin-flip symmetry. In contrast, the staggered magnetization operator is odd under such a transformation. Regarding the SYK Hamiltonian, spin-flip (i.e., 2i ↔ 2i − 1) turns one disorder realization into a different one with the same probability of occurring. As a consequence, given an instance of the SYK Hamiltonian that produces a certain dynamics ofR k , with k odd, there will always exist another realization that generates dynamics of equal amplitude, but opposite sign, leading to a vanishing ensemble average. N = 8, 12, 16, and 20 complex fermionic modes, and investigate the dynamics in the number preserved sector of half filling, for which we employ a state-of-theart, highly optimized exact diagonalization method. For the largest system size, the Hamiltonian matrix dimension is D = N !/[(N/2)!(N/2)!] = 184756. Due to the disorder, no symmetries other than particle number conservation can be used, and due to all-toall connectivity the matrix is denser than for models with finite-range interactions. However, thanks to the self-averaging nature of the model [63], with increasing N smaller numbers of disorder realizations suffice for satisfactory convergence [149] (we consider 90000, 2700, 400, and 50 ensemble members for the above values of N ). We note that, as an alternative to exact diagonalization, semiclassical methods have also been shown to successfully capture the dynamics of operator expectation values in Ref. [66].
With increasing N , faster equilibration as well as an approach to convergence is observed (see the right inset in Fig. 4, which highlights the dynamics in the transient time domain). A similar feature has been seen in the initial dynamics of other quantities for time evolution under the SYK [57] and other disordered, chaotic Hamiltonians [150,151], as well as random matrices [152]. The dependence on the system size can again be understood from the spectral analysis of the Liouvillian (see Sec. 5.3). For smaller N , the curves show oscillations before equilibrating to the steady state value. In addition, for N = 8, at large times the equilibrated curves slowly drift from the steady state value, with a rate that depends on the considered initial states. As a consequence, the rescaled curves cross zero and become negative at intermediate time, which is in accordance with Eq. (3) (see the left inset in Fig. 4, which highlights the approach of the curves to the steady state value). Both the transient oscillations and the drift can be attributed to finite-size effects, which become less prominent with increasing N . The oscillations may have random matrix origin, and are, in fact, also a feature of the SYK 2 model [111]. In contrast, the drift is due to finite-size induced non-exact uniformity in the distribution of initial states' amplitudes over the SYK 4 Hamiltonian spectrum. A similar drift has also been reported in the dynamics of the purity under Poissonian and Gaussian Unitary Ensemble random matrices [96]. The drift is reminiscent of the correlation hole, which appears in the evolution of survival probability, inverse participation ratio, and correlation functions [153] of chaotic systems. Alternatively, the slow drift constitutes the "ramp" of the dip-ramp found in the dynamics of the spectral form-factor [57] of SYK and random matrix models. Both terms describe the same phenomenon, which arises due to the long-range rigidity in the Hamiltonian spectrum. The depth of the correlation hole for equal time correlation functions is known to be suppressed for larger system size [153]. This, in another way, substantiates that the drift seen for N = 8 is a finite size effect for the observables considered in the present study.
To scrutinize the universality in more detail, we depict in Fig. 5a the evolution of G (E [M 2 ]) on logarithmic scales around the transient domain and at the verge of attaining the steady state (with N = 12 and 10000 disorder realizations). The curves show universality within the numerical precision, which until a time of about Jt ≈ 6 corresponds to the thickness of the curves. The vast part of the corresponding universal curve is well approximated by a Gaussian, which reaches well into intermediate times 1 Jt 10. This super-exponential behavior is followed by a marginal time domain with an indication of power-law decay to the steady-state. From Fig. 5a, one sees that the universality manifests throughout these regimes, while we cannot make predictions at larger times due to lacking convergence in disorder averaging. In Fig. 5b, we show the rescaled curve for N = 12, 16, and 20 to illustrate the system size dependence (choosing the initial FH ground state for U/J = 4 as a representative case). With increasing system size, the time-interval of the power-law decay seems to diminish, and the Gaussian appears to fit the curve for larger times (see inset). Closer to the steady state, a larger number of realizations is required for disorder averaging to converge. We further observe that the noninteracting FH ground state (U/J = 0) requires a larger sample size for convergence than other initial states.

Numerical results from master equation
In this section, we apply the open-system formalism of Sec. 4 to the SYK 4 Hamiltonian. We explicitly show the form of the jump operators and dissipation rates that were used for the ME simulations of the QFI presented in Fig. 2 and for the operator moments of the preceding section. To do so, one simply needs to rewrite the Hamiltonian of Eq. (1) in the generic form of Eq. (8), and then read off the jump operatorsĥ (α) lα and disorder functions ξ (α) lα that govern the dissipation rates as follows.
First, since Eq. (1) is a purely disordered Hamiltonian, we haveĤ 0 = 0 and the Liouvillian in Eq. (14) generates purely dissipative dynamics. Second, for the SYK 4 Hamiltonian, we have the multi-index l α = i 1 i 2 ; j 1 j 2 , and can identify three Hamiltonian termŝ H α with jump operatorŝ and corresponding time-independent disorder coeffi- Explicitly, for α = 1 the multi-indices are l 1 = i 1 i 2 ; i 1 i 2 with i 1 > i 2 , whereas for α = 2, 3 the multiindices are l α = i 1 i 2 ; j 1 j 2 with i 1 > i 2 , j 1 > j 2 and (i 1 , i 2 ) = (j 1 , j 2 ). Finally, we use the above expressions to determine the dissipation rates of Eq. (15). The relevant time integral is trivial in this case, and the rates are given by  The main point is that the time-independent disorder correlations of the SYK 4 model-or indeed any SYK q model given by Eq. (26)-yield dissipation rates in the Bourret-Markov ME that grow linearly in time. This property yields the super-exponential approach to equilibrium, as already discussed in Sec. 4 in the context of the QFI.  Figure 6 shows a comparison of ME and ED simulations for moment M 2 of operatorR defined in Eq. (6). As for the QFI, the super-exponential approach to equilibrium is captured, as well as the early and late time dynamics. We again observe a discrepancy between ED and ME simulations at intermediate times, a signature of correlation and memory effects as discussed in Sec. 4. As mentioned above, the ME curves contain no fit parameters. Recalling that hermitian jump operators guarantee the infinite temperature ensemble to be a (not necessarily unique) steady state of the general Liouvillian in Eq. (14), we determine the infinite temperature steady-state value of M 2 within the half-filling sector N = 2Q to be As Fig. 6 shows, this value agrees well with the steadystate plateaus of the exact unitary dynamics averaged over disorder realizations.

Analyses of Universality in ME formalism
Having used the ME framework to formally demonstrate the Gaussian decay via Eq. (17), we now study the origins of the observed universality within this formalism. To this end, we numerically obtain the spectrum and eigenmodes of the SYK 4 Liouvillian superoperator via a matrix representation thereof (see, e.g., Refs. [145,154,155] for detailed descriptions of such a procedure and for spectral properties of Liouvillian superoperators). As we will see, the population of various initial states and observables in the corresponding eigenspaces conspires to produce a universal curve under the rescaling G defined in Eq. (3).
In general, a superoperator L is not normal, and thus has distinct left and right eigenmodes. However, as a result ofĤ 0 = 0, (ĥ kα,lα ∈ R, the Liouvillian of the SYK model L(t) = 2tD is Hermitian and thus normal, as made explicit by the matrix representation given in Eq. (39). Therefore, the left and right eigenmodes of D coincide, and one can always form a Hermitian basis for each eigenspace. We use the index i ≥ 0 to label these eigenspaces, which in general have a d i -fold degeneracy. The d i Hermitian eigenmodes within the ith eigenspace are denoted asρ i,αi , where α i = 1, . . . , d i . The eigenmodes are orthogonal with respect to the Hilbert-Schmidt norm tr ρ † i,αiρ j,αj = δ i,j δ αi,αj and thus form a basis of B(H), the space of linear operators acting on H. For all (in our case typically degenerate) eigenspaces i, the corresponding eigenvalue λ i is real and negative. So, we order the eigenspaces according to the magnitude of their respective eigenvalues as |λ 0 | < |λ 1 | < . . . .
We can decompose any initial state and observable in B(H), respectively, aŝ Since Liouvillian dynamics are trace preserving, we have λ 0 = 0. Consequently, lim t→∞ρ (t) is given in terms of the eigenmodes corresponding to λ 0 [145]. For the Liouvillian of the SYK 4 model, we find λ 0 to be non-degenerate, implying a unique steady-state in the present case of study, in agreement with rather general conditions [156]. The Liouvillian spectrum {λ i } sets the time-scales of the dynamics of any observable quantity.
To analyze the universality of operator moments and the pure-state QFI observed under the rescaling G defined in Eq.
where In both cases, there exists only one non-zero amplitude A i * , and A i = 0 ∀i = i * , such that Eq. (25) reduces to the same Gaussian curve G Ô (t) = e −t 2 |λ i * | for all initial states.
We apply the above criteria to the first four moments M k of the staggered magnetizationR [see Eq. (6)] in the SYK 4 model. Table 1 lists the spectrum of D for a system of N = 8 modes at half filling. The corresponding distributions of c i,αi o i,αi and A i for M 2 are shown in Fig. 7 for different initial FH ground states. Since-for i > 0-c i,αi o i,αi = 0 only for i = 2, universality of type (i) for M 2 follows immediately with i * = 2. The reason why this occurs lies in the decomposition of the observable: We observe o i,αi = 0 for i > 0 and i = 2. This demonstrates why the choice of the initial state is completely irrelevantregardless of the values of c i,αi , there can be only a single A i = 0 for i > 0.
In Fig. 8, we display the behavior of moments M 1 , M 2 , M 3 , M 4 for the FH ground state at U/J = 10. For odd moments, we find A i = 0 ∀i, making them trivially universal: For M 1 , only c 1,α1 o 1,α1 = 0, whilst for M 3 additionally c 3,α3 o 3,α3 = 0. In either case, these terms are distributed near symmetrically about 0, such that the effective amplitudes vanish (for more details see Sec. 5.1).
In contrast, even moments have non-zero effective amplitudes in at least one eigenspace besides that of the steady state. In analogy to Fig. 8, we have verified for a range of initial FH ground states, as well as others such as the Neel state, that for a given even moment, any non-zero amplitudes A i always occupy the same eigenspaces. Concretely, for k ≥ 4, we find the same two non-zero effective amplitudes A 2 , A 3 , yielding an approximately universal super-exponential de- In summary, we find (i) odd moments vanish for all t (and are thus trivially universal), (ii) the second moment exhibits truly universal super-exponential decay as is shown in Fig. 7, and (iii) higher-order even moments display approximate universality. As this analysis shows, the Bourret-Markov ME reproduces the universal features observed in our exact numerics.
To conclude this section, we study the dependence of the Liouvillian spectrum on the system size N . We find that all non-zero eigenvalues decrease as N is increased from 6 to 8. In particular, for λ 2 -the only timescale entering M 2 (see Fig. 7)-we find, respectively, the values −0.28 and −0.33. This explains the shift to a faster equilibration time of M 2 with increasing N , as observed previously in Fig. 4. Note that this shift of M 2 exhibits a convergence, i.e., decreases as N is increased. We thus expect the eigenvalues λ i to not decrease indefinitely with N , but to individually approach some asymptotic value (for SYK q=2 this behavior can be shown analytically, but remains to be shown for q ≥ 4 [111]). However, the enlarged Hilbert-space, inherent to the process of mapping D to the matrix form of Eq. (39), limits our study of the Liouvillian spectrum to N ≤ 8, thus preventing us from probing this convergence.
In this context, we emphasize that the aim of the proposed ME framework is not to provide an efficient way of simulating the disorder-averaged dynamics. Rather, the aim is to gain additional insights by establishing a theoretical mapping to an open quan-tum system. For example, our present study reveals a non-trivial highly-degenerate eigenspace structure of the effective Liouvillian superoperator, through which we can explain the observed universality. The origin of this eigenspace structure is, briefly put, an operator size symmetry where each degenerate eigensector corresponds to a set of operators sharing a common operator size. For an in-depth study of this structure, we refer the interested reader to our follow-up work Ref. [135].   Fig. 7, but for different moments M1, M2, M3, M4 ofR, for the initial FH ground state at U/J = 10. For the highest moment, occupation of two eigenspaces can be observed, indicating that universality deteriorates in many-body operators.

Conclusion and Discussion
In summary, we have theoretically investigated post quench equilibration dynamics of a system of randomly interacting fermions described by the complex SYK 4 model. By numerically studying the disorderaveraged exact evolution of a set of local observables and their higher-order moments, we find that the equilibration process is universal. The curves illustrating the equilibration of different initial states overlap throughout the dynamics under a straightforward rescaling, revealing the independence of the dynamics on chosen initial states. The equilibrated steady state, which is the Gibbs infinite temperature state in the present study, is reached in fast time-scales of leading-order processes determined by the variance of the disordered interaction. In addition, the universal equilibration curve can be well approximated by a Gaussian, yielding fast super-exponential equilibration dynamics.
In order to achieve an analytical understanding of the numerical findings, we have formulated a theoretical framework based on the Novikov-Furutsu theorem. This framework describes how a disordered quantum many-body system undergoes an effective dissipative dynamics due to phase mixing in the ensemble averaged evolution rather than to interactions with a heat bath [96,140]. Thanks to the generality of the formulation, its scope for applications extends beyond the present investigation. Employing Bourret and Markov approximations, we obtain a Lindblad master equation that successfully describes the key features of the observed super-exponential equilibration dynamics under the SYK 4 model. Furthermore, a spectral analysis of the corresponding Liouvillian superoperator illuminates the exact universality of loworder moments, representing few-body observables, as well as an approximate universality of higher-order moments representing many-body observables.
The Novikov-Furutsu theorem has been used extensively in the literature for the study of systems with noise of short correlation time [97][98][99], and equations equivalent to those derived in Sec. 4 have been obtained to second order in perturbative noise strength [100-102, 138, 139]. In the present scenario, where noise correlation times are formally infinite and where the disordered interactions provide the dominant (because only) energy scale, there is at first sight no reason for such perturbative approaches to hold. Yet, the strong chaoticity of the SYK 4 model leads to a fast decorrelation, making the Bourret-Markov approximation an excellent description of the exact dynamics.
The salient features of the universal curve occur on sizable absolute scales and very fast time scales, on the order of the mean interaction strength, and they can be extracted from the observation of local observables following a simple global quench. Thus, the discussed effects should be readily observable in forthcoming laboratory implementations of the SYK model, for which several proposals have recently been put forward [71][72][73][74][75][76]157]. Moreover, the mapping to the purely dissipative Lindblad equation may also open new ways for simulating SYK matter using engineered open quantum dynamics.
We hope our findings will also stimulate further theoretical investigations to obtain and understand universality in out-of-equilibrium dynamics of other quantum many-body systems. A topic of our immediate future investigation is, e.g., in how far the universality survives in disordered models without allto-all connectivity. Further, while we have considered in this work only the SYK 4 model, our findings hold equally well for other SYK q models [111]. In the future, it will also be interesting to apply our master equation approach to other disordered models. Whilst time-independent Gaussian disorder will always yield a Gaussian decay according to our formalism, the universality (initial state independence) would depend on details of the Hamiltonian. In particular, in models with a clean contribution,Ĥ 0 = 0, where one can expect a complex interplay between disorder and clean dynamics, such as the many-body localization transition [158][159][160][161]. Various mathematical extensions of the presented master equation framework are also possible, e.g., to treat non-Gaussian disorder [134,148]. Finally, as this framework can naturally include dephasing noise with arbitrary correlation spectrum, it permits one to estimate the interplay between dissipation due to an external environment and dephasing due to disorder averaging, which is central, e.g., for environmental assisted quantum transport [162][163][164][165].

Acknowledgments
We gratefully acknowledge useful discussions with Jean-Philippe Brantut, Julian Sonner, and Ricardo Costa de Almeida. We acknowledge CINECA HPC project ISCRA Class C: ISSYK-HP10CE3PVN and ISSYK-2 (HP10CP8XXF). We acknowledge support by Provincia Autonoma di Trento and the Google Research Scholar Award ProGauge. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 804305). This project has benefited from Q@TN, the joint lab between University of Trento, FBK-Fondazione Bruno Kessler, INFN-National Institute for Nuclear Physics and CNR-National Research Council.

A Details of SYK q model
Here, we present the Hamiltonian of the general SYK q model of which Eq. (1) is a special case with q = 4. The SYK q Hamiltonian is governed by disordered allto-all (q/2)-body interactions, and reads [48,112] where K = [(q/2)!(q/2 − 1)!]/N q−1 ensures the extensivity of the model. The interaction strengths in Eq. (26) are complex Gaussian random variables, i.e., the real and imaginary parts of {J i1...i q/2 ;j1...j q/2 } are independent and normally distributed with variances parameterized by J ∈ R >0 as Furthermore, the interaction strengths satisfy where P and P perform permutations of the indices, and sgn(P), sgn(P ) = ±1 denote the sign of the permutations. The first equality ensures Hermiticity of the SYK q Hamiltonian, whereas the second is due to the fermionic anticommutation relations of the creation and annihilation operators in Eq. (26).

B Higher-order moments
Here, we present additional results for the quench dynamics of moments k = 8, 10, 12 (see Fig. 9) of the staggered magnetizationR defined in Eq. (6). The rescaled dynamics (right column) show approximately universal Gaussian equilibration similar to the other moments k < 8 studied in the main text. However, unlike the case of k = 2, the universality is not exact for these higher-order moments, as is explained in detail in Sec. 5.3. In addition, the rescaled curves shift towards earlier time with increasing k, which is evident from the lower insets in Fig. 10. Further, the rescaled curves seem to converge at sufficiently large k. In order to explain this characteristic, we plot the ratio A 3 /A 2 of the effective amplitudes [defined below Eq. (25)] between the occupied eigensectors of the Liouvillian (see top inset of Fig. 10). The ratio shows a monotonic increase and a saturation with respect to k.
Together with the fact that |λ 3 | > |λ 2 | (see Table 1), this explains-in accordance with Eq. (25)-the shift of the rescaled curve to earlier times, as well as its convergence, with moment order. We note that this convergence with moment order, and its explanation in terms of effective amplitude ratios, holds also for other choices of observables and/or initial states (not reported here), though the direction of convergence need not always be to earlier times but can also be to later times. Finally, we comment on the lack of exact universality observed for moments with k > 2: The kth moment ofR contains interactions of up to k fermionic modes, and thus probes increasingly nonlocal physics for larger values of k. As a limiting case of a truly global many-body quantity, we show in Fig. 11 the disorder-averaged survival probability (or fidelity) E | ψ(0)|ψ(t) | 2 . For this quantity, the curves corresponding to different initial states do not collapse. Whilst this suggests that universality is absent for (superpositions of) highly non-local observables, we do not exclude the possibility that one may be able to specifically construct such an observable which does exhibit universality. Finally, we note that the dynamics of the survival probability, namely an initial Gaussian decay followed by oscillations at intermediate times, is in accordance with the well established behavior of the post-quench dynamics of the survival probability within random matrix theory [130,150].

C Additional observables
Here, we extend our investigations to other observables in order to show the generality of our key findings. For this, we first consider the following 4-local operators defined in terms of theκ i s introduced in Eq. (6), (29) The system average of these operators is defined aŝ jŜ j , where j ∈ {1, 2, ..., N/4}. For an N = 16 system, we can construct four such 4-local operators, i.e.,Ŝ 1 ,Ŝ 2 ,Ŝ 3 , andŜ 4 . In Fig. 12, we show the representative evolution of S 1 (t) = ψ(t)|Ŝ 1 |ψ(t) as well as of S(t) = ψ(t)|Ŝ|ψ(t) . As in the cases of the QFI and moments of the staggered magnetization, we recover the super-exponential universal equilibration dynamics of these 4-local operators.
All the operators considered so far are diagonal with respect to the Fock-space spanned by the occupation number basis vectors {|n 1 , n 2 , ..., n N }. Here, we additionally investigate the dynamics of a non-diagonal For all curves the initial state is U/J = 10, and the dark to light shading corresponds to different moment orders k = 2, 4, 6, 8, 10 and 12, respectively. A shift to earlier times with increasing k is evident, and is further highlighted by the two lower insets. There is also an indication of convergence to a fastest decay time with increasing k. This is further studied in the top inset, which shows the ratio A3/A2 of the effective amplitudes of the kth moments in accordance with Eq. (25). The monotonic increase, and indication of saturation, with increasing k shows that the faster time-scale |λ3| is favored with increasing moment k.  Figure 11: Ensemble averaged survival probability for system sizes N = 8 (green), 12 (yellow), 16 (blue), and sample sizes 90000, 10000, 400, respectively. Shadings of a given color represent different initial states, as in Fig. 9. There is a clear lack of universality with respect to different initial states from early to late times. The early-time Gaussian decay followed by oscillations at intermediate times are expected from random matrix theory (RMT) [130,150]. In contrast to RMT where universality is expected for the survival probability, here we observe an initial state dependence, which can be attributed to the structure of the SYK Hamiltonian, imposed by the fermionic statistics.
In Fig. 13, we present the dynamics of the QFI, F Q , computed with respect to this operator. Similar to all the previous cases, we again obtain the superexponential universal equilibration dynamics of this observable. shading represents different initial states, as in Fig. 9. Again, the rescaled data (right column) is well fitted by a Gaussian, exp −(Jt/τ ) 2 , with τ = 1.52.

D Details of ME derivation
Here, we provide further details on the derivation of the Lindblad master equation (ME) presented in as δρ [ξ, t] δξ (α) in which the step-function Θ arises from causality. The above recursive expression for the functional derivative yields a series of nested commutators, which to lowest order reduces to the first term of Eq. (32). As mentioned in the main text, truncation to this lowest order is motivated by the fact that the resulting evolution equation [given by Eq. (13)] is formally equivalent to that obtained when making the decorrelation assumption in the study of stochastic evolution equations [136,137]. This is readily seen in the interaction picture generated byĤ 0 (t): Integrating the von Neumann equation of an individual Hamiltonian realization one obtains a selfconsistent integral equation forρ(t). By inserting this back into itself once, taking the disorder average and finally differentiating with respect to time, one finds an equation which reduces to the interaction picture version of Eq. (13) after performing the decorrelation approximation E ξ (α) Lindblad form.-The final Lindblad master equation of Eq. (14) is obtained from Eq. (13) by first making the Markov approximationρ(t ) ρ(t) and then expanding the double commutator. Simplifying our notation, this expansion is l,k f l,k (t)(ĥ lĥkρ −ĥ lρĥk −ĥ kρĥl +ρĥ kĥl ), (33) which is already reminiscent of a master equation in standard form. The latter is obtained in a final step in which we require the correlations to be symmetric under an exchange of the indices, i.e., f l,k (t) = f k,l (t). This is trivially fulfilled for static processes such as those of the SYK model. For continuous processes, our requirement is equivalent to symmetry in time, . We can then regroup the terms of Eq. (33) as We thus finally obtain the Lindblad master equation in non-diagonal form, given in the main text by Eqs. (14)- (15).

Cross-correlations of dissipation rates.-
Here we comment on the origin of the cross-correlations for α = 2, 3 which exist in the dissipation rates, defined in Eq. (21), of the SYK 4 model. To rewriteĤ SYK4 in the form of Eq. (8), we partially order the indices as i 1 > i 2 , j 1 > j 2 and use the anti-symmetry [see Eq. (28)] of the SYK interactions. We then find the three disorder terms given in Eq. (20). The point is now that for α = 2, 3, correlations exist between pairs ξ kα if multi-indices l α , k α are mirror images, i.e., if l α = i 1 i 2 ; j 1 j 2 and k α = j 1 j 2 ; i 1 i 2 . This is because ReJ i1i2;j1j2 = ReJ j1j2;i1i2 and similarly ImJ i1i2;j1j2 = −ImJ j1j2;i1i2 , due to Eq. (28).
Limit of time-independent processes.-Here we explicitly show that the Lindblad ME of Eqs. (14)- (16), formally derived for Gaussian processes with arbitrary time-dependence, is valid in the limit of timeindependent (static) processes ξ We now have modified correlations E ξ (α) lαρ (t) , in whichρ(t) is a function (as opposed to a functional) of the Gaussian distributed random numbers ξ (α) lα (as opposed to random functions). These correlations may, however, still be treated via the Novikov-Furutsu theorem, which simplifies accordingly to [134] E ξ where now we have disorder correlations F (α) with infinite correlation time. This static version of the Novikov-Furutsu theorem is also known as Stein's lemma [166].
The right-hand-side of Eq. (36) now contains an ordinary total derivative (as opposed to a functional derivative), which we again obtain to lowest order from the integrated von Neumann equation as Substituting Eqs.
This is exactly the evolution equation that one obtains by taking the limit of time-independent processes, i.e., of infinite correlation times in Eq. (13). lα,kα . This shows that starting from time-dependent noise and then taking the limit of infinite correlation time is equivalent to working with static noise all the way. The advantage of the time-dependent formulation used in the main text is that it naturally incorporates static disorder and temporally fluctuating noise on the same footing.

Matrix
representation of Liouvillian superoperator.-We numerically obtain the spectrum {λ i } and eigenmodes {ρ i } of a Liouvillian superoperator L by numerically diagonalizing a matrix representation thereof. This representation is obtained by mapping L →L, withL an operator acting on the duplicated Hilbert space H ⊗ H (see for instance Ref. [145]). Under this map, a density matrixρ becomes a vector | ρ ∈ H ⊗ H obtained by stacking the columns of the matrix representation ofρ. Left and right multiplication with operators transform asÂρB →B ⊗Â | ρ , whereB denotes the transpose ofB. For our general Liouvillian [see Eq. (14)], we then have the matrix representation We obtain the matrix representation of the SYK 4 LiouvillianL(t) = 2tD (factoring out 2t as in the main text) by settingĤ 0 = 0 and inserting Eqs.

E Details on numeric simulation of exact dynamics
Here, we provide a brief description of the algorithm implemented to solve the exact dynamics reported in this manuscript. We exploit the particle number conservation of the SYK HamiltonianĤ SYK4 [see Eq. (1)] by restricting our simulations to a given particle number sector of the Hilbert space. In particular, all simulations reported in this manuscript are performed within the half-filling sector where the fermion number is Q = N/2, and the Hilbert space dimension is D = N !/((N/2)!) 2 (see Sec. 2).
The matrix representation ofĤ SYK4 is constructed with respect to the fermion mode occupation number Fock basis of the Q = N/2 sector: Each Fock state |s a , where a = 1, . . . , D, is represented by one of the N -bit strings s a of which Q bits are 1 in order to represent the occupied fermion modes. The Hamiltonian's matrix elements H ab ≡ s a |Ĥ SYK4 |s b are non-zero only for those pairs of states whose bit string representation have a Hamming distance d(s a , s b ) = 0, 2, 4. This follows from the quartic operatorsĉ † i1ĉ † i2ĉ j1ĉj2 appearing inĤ SYK4 , where d(s a , s b ) = 0 corresponds to (i 1 , i 2 ) = (j 1 , j 2 ), d(s a , s b ) = 2 to i k = j l , i k = j l for k, k , l, l = 1, 2, and d(s a , s b ) = 4 to i 1 = i 2 = j 1 = j 2 . These non-zero elements are populated by independent random complex Gaussian variables J i1i2;j1j2 , in accordance with Eqs. (27) and (28) to ensure Hermiticity of the Hamiltonian and antisymmetry of the interaction amplitudes under permutation of the indices. In this way, a single realization of the random SYK 4 Hamiltonian is constructed.
For system sizes N ≤ 14, D is sufficiently small such that the full eigensystem (energy basis) of the above matrix representation can be solved exactly. For this we employ diagonalization routines from LAPACK [167]. Time evolution of a given initial state |ψ(0) is then solved by rotating the corresponding vector representation into the SYK energy basis and calculating the time-dependent phases exp(−i n t) for any time t ∈ R, where n for n = 1, . . . , D are the eigenenergies ofĤ SYK4 .
For system sizes N > 14, D is so large as to prohibit the above exact numeric solution of the entire eigensystem. Instead, we utilize a Runge-Kutta 4 (RK4) method to solve the Schrödinger equation forĤ SYK4 . The matrix-vector multiplication employed within the RK4 method makes use of a sparse matrix representation ofĤ SYK4 in order to exploit the large amount of zero matrix elements, and thus enhance the speed of the algorithm. Further reduction of computation time is achieved by parallelizing this matrix-vector multiplication via MPI methods [168], allowing us to exactly solve (within numeric precision) the dynamics for system sizes up to N = 20 at Q = N/2. To ensure accuracy of this RK4 based algorithm, we benchmark its dynamics against those of the above exact diagonalization scheme for N ≤ 14, verifying that they agree.
Finally, we average over the dynamics of multiple disorder realizations ofĤ SYK4 . To this end, we again utilize MPI methods to solve the dynamics of multiple disorder realizations in parallel.