Effective dynamics of disordered quantum systems

We derive general evolution equations describing the ensemble-average quantum dynamics generated by disordered Hamiltonians. The disorder average affects the coherence of the evolution and can be accounted for by suitably tailored effective coupling agents and associated rates which encode the specific statistical properties of the Hamiltonian's eigenvectors and eigenvalues, respectively. Spectral disorder and isotropically disordered eigenvector distributions are considered as paradigmatic test cases.


I. INTRODUCTION
Disorder is the expression of a lack of knowledgee.g., on a physical system's conformation or of a potential landscape, on scales that affect the system evolution. The characterization of a disordered system therefore requires a statistical approach, since reliable predictions on reproducible features of the system behavior can only be extracted by averaging over suitably chosen distributions of those uncontrolled properties. The specific choice of the distribution to be averaged over has potentially strong impact on the predicted behavior. In turn, the observed behavior may provide strong, though not necessarily unambiguous, hints on the underlying disorder's structure, as well known, e.g., from the classical Hamiltonian flow in mixed phase spaces (where one seeks a statistical characterization of the time evolution of ensembles of initial conditions) [1][2][3].
On the quantum level, the disorder average has yet another important consequence: It implies an average over the accumulated phases associated with the eigenstates of every realization of the underlying random Hamiltonian (see Fig. 1). In general, this induces a loss of phase information, hence decoherence, while, simultaneously, dramatic disorder-induced interference effects may prevail under the disorder average. Note that, here, decoherence is an immediate consequence of the disorder average, and must not be mistaken for an irremediable loss of information to a large environment (bath) as in the context of open system. The "lost" information is encoded in higher-order correlations such as intensity correlations of a speckle potential [4], or in fluctuations of the measured conductance across individual samples [5]. Such intricate indicators of coherence effects may, however, be difficult or even impossible to measure if single realizations of the disordered potential cannot be reliably resolved for several successive experimental runs as needed for the measurement of an observable. As an illustration, think of the propagation of photons in the presence of a turbulent atmosphere [6]. It therefore arises as a natural question * Chahan.Kropf@physik.uni-freiburg.de † A.Buchleitner@physik.uni-freiburg.de what we can learn on the underlying disorder from the time-evolution of the ensemble averaged state alone, and especially from the coherent and incoherent content of the latter.
Prominent examples of robust interference effects which prevail even under the ensemble average are localization phenomena in the quantum transport theory of disordered or chaotic, dynamical systems [7][8][9][10]. So far, however, these are discussed in terms of the associated spectral properties or of the signature of quantum interference in characteristic, experimentally observable quantities -contrasted against classical predictions. Scattering theoretical approaches account for the above phase average and identify robust interference contributions, though almost always in the stationary state. Furthermore, they only rarely quantify the relative weight of those coherent contributions which survive as compared to those which are eliminated under the disorder average, let alone the dynamical evolution of this ratio until stationarity is reached [11,12].
No formalism is so far available which allows to assess the effective dynamical quantum evolution of an arbitrary initial state in real time, and to directly associate characteristic time scales and couplings thereof with the underlying disorder. Note, however, that this represents a natural and -given the wealth of disordered quantum transport problems -substantial expansion of the generic playground of the theory of open quantum systems, by substituting an uncontrolled environment by a static, operator-valued random perturbation of the system Hamiltonian. The relevance of such generalization of open system theory appears highly plausible, since it will pave the way for a systematic assessment of the question which type of disorder allows to exploit which type of quantum coherent phenomena, on transient and/or asymptotic time scales. Indeed, the direct experimental monitoring of the dynamical evolution of disordered quantum systems has now moved into reach [13][14][15][16], and it therefore appears timely to develop tools to distinguish coherent and incoherent [11,12] features thereof.
Here, we derive effective dynamical evolution equations for the ensemble-averaged state of disordered quantum systems, in the form of master equations, which, by their very structure, precisely meet the above purpose, and The ensemble averaged state ρ(t) is then obtained by averaging over all ρi(t). The different realizations of the generating Hamiltonian may, e.g., result from slow parameter drifts between subsequent runs of an experiment, or from microscopically distinct potential landscapes in macroscopically identically prepared experimental settings.
show how the statistics of the disorder enters the unitary part as well as the Lindblad operators and associated rates, as the equation's specific ingredients. We apply this theory to exemplary cases of quantum systems with random spectra and randomly distributed eigenvectors for which our method can be applied without any approximations.
The article is structured as follows: In Section II, we introduce disordered quantum systems and the dynamics which emerge from the ensemble average. In Section III, a direct derivation of a master equation valid in the limit of short times is presented (see also [17]). Section IV presents a more general method to obtain disorder master equations at all times. We then elaborate these methods for a single qubit with spectral disorder in Section V, investigating the impact of different eigenvalue distributions. We continue in Section VI with the generalization of spectral disorder to higher dimensions and study the role of correlations among eigenvalues. Finally, we describe in Section VII the resulting master equations for unitarily invariant disorder ensembles, with which we illustrate the effect of disorder in the eigenvectors, and evaluate them for Gaussian and Poissonian eigenvalue statistics. Section VIII concludes the paper.

II. DISORDER ENSEMBLE AVERAGE
To start with, let us briefly introduce quantum systems with disorder on the level of the Hamiltonian, and the corresponding ensemble average dynamics. Basic properties of the latter are exposed and their description in terms of quantum master equations outlined.

A. Single realizations and ensemble average
We consider an isolated, disordered quantum system of dimension d. The disorder may be characterized by an ensemble of time-independent Hamiltonians H λ , occurring with probability p λ , where the (multi-)index λ labels the different realizations.
For each realization H λ of the disorder, the corresponding state ρ λ (t) follows the von Neumann equation of motion, whereρ λ (t) = (d/dt)(ρ λ (t)). Its solution is given by The initial (t = 0) state ρ(0) is taken to be identical for all realizations.
The ensemble average state ρ(t) (all ensemble average quantities will be marked with a bar) is obtained by the weighted sum over all realizations (for convenience denoted by an integral throughout), For a given observable B = b b |b b|, the ensemble average state delivers the average probability of measurement outcome b, which is given by p(b) = dλp λ Tr[|b b|ρ λ ], directly as Formally, the dynamics of the ensemble average state (3) are described by a family of linear maps Λ t , parametrized by the time t, from the set of density matrices D onto itself: In the following, we will describe some characteristic properties of this family of maps.

B. Properties of the dynamical map
The dynamics of each single realization is unitary and thus Hermiticity-preserving, trace-preserving, and completely positive. The linearity of the averaging procedure (3) then implies that Λ t also has the latter three properties. Hence, the ensemble average dynamics describe by construction legitimate quantum dynamics [18]. But, in contrast to single realizations, the ensemble average dynamics are in general nonunitary (as we will explicitly show below), implying that Eq. (3) cannot be subsumed by a single unitary operator.
A general property of the dynamical map is that the maximally mixed state is an invariant, Λ t [1] = 1, indicating that the map is unital or bistochastic, in mathematical and statistical physics). In that sense the dynamics emerging from ensemble averages are more restricted than those of open systems, which can in general also be non-unital, as, e.g., in spontaneous decay processes of atoms.
As a direct consequence of the non-unitarity, purity, Tr ρ 2 , is in general not conserved. More precisely, the unital map Λ t [ρ(0)] is always purity-decreasing, Tr ρ 2 (t) ≤ Tr ρ 2 (0) [19]. As we will see below, this does however not imply that purity decays monotonously; it can increase locally, although it can never exceed the purity of the initial state.

C. Master equation description
For formal clarity, we here summarize several key properties of quantum master equations in general [20], and especially in the context of disorder dynamics.
As was briefly outlined in the introduction, we wish to describe the dynamics of the disorder ensemble average, formally defined by the map (5), in terms of a master equation. Since the ensemble average dynamics in general exhibit decoherence, they cannot be captured solely by the von Neumann equation of motion (2) with some appropriately chosen Hamiltonian. Quantum master equations generalize the latter towards non-unitary dynamics while still guaranteeing the preservation of the Hermiticity and the trace of the state, and maintaining the complete positivity and the linearity of the underlying dynamical map. The general structure of quantum master equations imposed by these properties is given by the Lindblad form [21,22] This form complements the coherent dynamics of the von Neumann commutator (2) by the incoherent dynamics induced by the Lindblad operators L k (t) and their corresponding decoherence rates γ k (t). The curly brackets {A, B} = AB + BA denote the anti-commutator. We remark that the Lindblad form (6) does not uniquely fix the Hamiltonian and the set of Lindblad operators [20]. Let us stress that any integro-differential master equation including a memory kernel can be cast into the time-local form (6) [23] (see also Appendix D). Possible memoryeffects in (6) are encoded in the time-dependence of the rates and Lindblad operators. We focus on this time-local form because it offers the physically most transparent interpretation for our purposes.
Our goal is to determine the relations between the composition of the disorder ensemble (1) and the corresponding emerging master equation (6), where the Lindblad operators and their rates will capture the collective (incoherent) dynamical impact of the disorder distribution. As we will show below, the ensemble average generically gives rise to a time-dependent Lindblad term, including negative decoherence rates (which are often considered as a signature of memory effects in the open system context [24]). While the latter can in general give rise to unphysical dynamics, in the considered case the emerging master equation will, by its very construction from the ensemble average, always be physically consistent. In the case of time-independent Lindblad operators and time-independent, positive rates, one usually speaks of a Markovian Lindblad master equation in the strict sense, describing a dynamical semi-group.

III. SHORT-TIME MASTER EQUATION
For short times one can derive a general expression for the ensemble average master equation. To see this, let us consider the time evolution of a single realization over a time step dt, ρ λ (dt) = U λ (dt)ρ(0)U † λ (dt), and expand the time evolution operators to second order. One then obtains As we will see, the expansion to second order is necessary, since the leading incoherent dynamical contributions only appear at this order. If we now take the ensemble average of (7), isolating the contribution of the average Hamiltonian H from the second order term, we obtain where we have introduced the Hermitian operators L λ = H λ − H. Eq. (8) can be conceived as the second-orderin-time expansion of a master equation in Lindblad form (6), which approximates the dynamics at short times (the step from Eq. (8) to Eq. (9) is not completely trivialfor more details see [17]): Each realization λ gives in general rise to a Lindblad operator L λ and a corresponding decoherence rate γ λ (t), where ω 0 is introduced as a characteristic energy scale of the system. We thus find that the ensemble average dynamics are described, at short times, by timeindependent Lindblad operators given by the deviations of the single realization Hamiltonians from the average Hamiltonian, and by linearly-in-time increasing decoherence rates. Note that the above derivation is based on the second-order in time expansion (7) and does not involve any further approximation, e.g., on the level of the ensemble averaging. As we will see below, the particular structure of the disorder (e.g., when garnished by symmetries) will often allow us to transform (9) into a significantly simplified master equation. Moreover, we will see that, in our cases, the short-time master equation already reflects the structure of the master equation at arbitrary times, with only the time dependence of the decoherence rates modified, while in other cases the finite-time dynamics may also give rise to modified Lindblad operators. The coherent part of the short-time dynamics is induced by the average Hamiltonian H. Note that, since the decoherence rate vanishes at t = 0, the Lindblad term (the second line in Eq. (9)) does not contribute to first order in time. Furthermore, we remark that in deriving the short-time master equation (9) we had to assume that the average Hamiltonian H exists; if this is not the case, the shorttime dynamics may differ from (9,10), as in the examples which will be discussed in Sections V A and V D below.

IV. FINITE-TIME MASTER EQUATION
We are not aware of a general expression, similar to the short-time master equation (9), for the ensemble average master equation at arbitrary times. In this section, we outline a method that can be employed in order to determine the master equation for arbitrary disorder distributions given the following holds: Let us assume that the inverse Λ −1 t [ρ] of the dynamical map Λ t [ρ] defined in Eqs. (3) and (5) exists. We can then write . A time-local differential equation is formally obtained by taking the time derivative of Eq. (5), where • denotes map composition. Note that this can be done for any sufficiently well-behaved dynamical map Λ t [ρ]. The issue of the possible non-existence of the inverse Λ −1 t [ρ] will be discussed at the end of this Section. In the following we describe an explicit method to compute the mapsΛ t and Λ −1 t , and how to obtain from this a master equation in Lindblad form. The method is based upon a matrix approach presented by Andersson and Hall [25,26]. It can in principle be applied to any kind of Hamiltonian disorder, but the computational complexity can be substantial, making approximations often unavoidable. Note, however, that for all the examples presented in this paper (see Sections V-VII), we can derive the disorder master equation without any approximations.
Expressed in terms of an orthonormal basis {|j } d j=1 , the ensemble average state (3) reads The basis {|j } d j=1 can be chosen suitably to ease the computation of the inverse of the average dynamical matrix defined in Eq. (13) below. For convenience we adopt a vector notation and define the d 2 × d 2 average dynamical matrix F (t) and the d 2 × 1 average density vector ρ component-wise by and ρ jk (t) ≡ j|ρ(t)|k , where (jk) and (rs) are double indices with j, k, r, s ∈ {1, 2, ...d}. We remark that the average dynamical matrix F (t) contains the same information as the d 2 × d 2 Choi matrix [27]. (Equivalent dynamical matrices with different index orderings exist in the literature [28,29].) In terms of the average dynamical matrix F (t) and the average density vector ρ, the ensemble average state (3) is obtained by the standard matrix product ρ(t) = F (t) · ρ(0). Based on this representation, the inverse and the time derivative of the dynamical map can be computed using standard matrix operations. Concretely, the differential equation (11) can now be written aṡ where the d 2 × d 2 matrix Q(t) ≡Ḟ (t) · F −1 (t) represents the mapΛ t • Λ −1 t . This implies, in terms of the components of Q and ρ: The final Lindblad form (6) is then obtained by expanding the operators |j r| and |s k| in a Hermitian operator basis and collecting the different terms using the hermicity of ρ(t) and trace-preservation, Tr ρ (t) = 0 (this step is also performed in the textbook derivation of the Lindblad master equation in open quantum systems theory [20]) . More precisely, one chooses a Hermitian operator basis {A m } By setting n = 0, the orthogonality relation implies Tr[A m ]| m =0 = 0. Any operator basis satsifying (16) can be chosen (note that the Hermicity condition is for convenience only). However, a suitable choice may simplify subsequent calculations. A natural choice are, e.g., the d 2 − 1 Gell-Mann matrices [30], which are a direct extension of the Pauli matrices to higher dimensions. The Gell-Mann matrices (here denoted with G in order to avoid confusion with the unspecified basis operators A) can be separated into sets of symmetrical {G s }, antisymmetrical {G a }, and diagonal {G d } matrices. Given the basis {|j } d j=1 , we have, with the normalization required by (16), Coming back to the general operator basis (16), Eq. (15) rewritten in terms of the A m yields (here and in the following we drop the time-dependence of ρ in the notation of the master equations) where we introduced for clarity, and as an intermediary step, the Hermitian matrix with the abbreviation A m,rj = r|A m |j , the rj component of A m . The terms with m = 0 or n = 0 are separately collected iñ From Eq. (18) and since Tr ρ = 0, one finds that If we now introduce the effective Hamiltonian we arrive at the non-diagonal Lindblad forṁ with the decoherence matrix Γ mn (t) = C mn (t), (m, n > 0, i.e., the matrix (19) without the first line and the first column), and the (Hermitian, time-independent) Lindblad operators A m . Note that the ensemble average (3), even when taken over the unitary dynamics arising from static Hamiltonians as considered in this article, gives rise to a possibly time-dependent effective Hamiltonian H(t), which does in general not coincide with the average Hamiltonian H. We emphasize that this time-dependence does not mediate the effect of an external driving potential as would be usual in the context of open quantum systems, but is a direct consequence of the composition of the disorder ensemble (1). An example is discussed in detail in Section V D.
Since the decoherence matrix is Hermitian, and γ k (t) the kth eigenvector and eigenvalue of Γ(t), respectively. This step requires explicit knowledge of the decoherence matrix Γ(t). We then obtain the diagonal Lindblad form (6) of the disorder master equation with Lindblad operators L k (t) = d 2 −1 m=1 V mk (t)A m and strictly real decay rates γ k (t). In general, the diagonalization leads to time-dependent and non-Hermitian Lindblad operators. However, below we will also give an example of a disorder distribution giving rise to a time-independent master equation.
Note that, in the spirit of the short-time expansion in Sec. III, one can expand the decoherence matrix Γ(t) in time, Γ mn (t) = Γ mn (0)+Γ mn (0)t+Γ mn (0)t 2 /2+. . . This then reexpresses the nondiagonal Lindblad master equation, where each order of the time expansion of Γ(t) leads to a Lindblad term, i.e., each time order can be diagonalized independently. When all derivatives d k dt k Γ mn (0) commute, there exists a set of time-independent Lindblad operators specifying the diagonal form (6) at all times. If one truncates the expansion after the linear term, one recovers the incoherent part of the short-time master equation (9) (provided it exists, i.e., provided H is well defined, which is not the case for certain disorder distributions, such as for example the Cauchy-Lorentz diagonal disorder considered in V A), however in a different representation.
The presented method to convert a dynamical map into a time-local master equation requires the existence of the inverse map Λ −1 t [ρ], cf. Eq. (11). A necessary condition for the inverse not to exist is that two different initial states evolve into the same state at a given finite time t > 0 [26]. As we will see below (e.g. in Section V), such coincidence can indeed occur for the ensemble average dynamics. However, this does not necessarily imply that one cannot find the corresponding master equation. If the inverse does not exist only at isolated points in time, one can still formulate a master equation which then reflects the non-existence of the inverse by diverging decoherence rates. In Section V, diverging rates will for example arise for a single qubit subject to uniform spectral disorder. As we will see, there the divergences are a consequence of the compact support of the disorder distribution.
Indeed, in all situations considered in this article, the inverse map exhibits at most isolated divergences. We conjecture that this is a general property of dynamical ensemble averages on finite-dimensional Hilbert spaces as defined in (3), due to the quasiperiodicity of the timeevolved state inherited from the discrete spectrum of U λ .
Finally, we emphasize that the method to derive a master equation starting from a quantum dynamical map Λ t [ρ] presented in this section is not restricted to ensemble average dynamics, but can be applied to any invertible dynamical map. Note however that, in some cases, other, more direct, methods, may be preferable. For instance, in the case of unitarily invariant disorder (cf. Section VII), or when considering random unitary channels as they are usually defined in quantum information theory (ρ = j p j (t)U j ρU † j ) [31], one can directly invert the map Λ t . The Choi-Jamilkowsky isomorphism [32,33] may provide us with yet another alternative method.

V. SPECTRAL DISORDER: SINGLE QUBIT
The first conceptual benchmark situation we focus on has the disorder occurring in the spectrum of the Hamiltonians, while all realizations share the same eigenvectors (for the moment we consider the general case of a d-dimensional quantum system), The disorder is then fully characterized by the distribution p λ of the eigenenergies ω 0 λ ≡ ω 0 (λ 1 , ..., λ d ) T , where ω 0 ∈ R + is the characteristic Larmor precession frequency. Such scenario for example describes an ensemble of non-interacting spins 1/2 in a static magnetic field with spatial inhomogeneities, or which fluctuates in intensity from one measurement to another, a situation for instance encountered in magnetic resonance spectroscopy [34] or in ion trap experiments [35]. We start by considering the case of a spectrally disordered qubit (i.e. a system of dimension d = 2), illustrating in some detail the matrix approach to compute the disorder quantum master equation outlined in Section IV. As we show, even for such a simple system, the ensemble average can exhibit nontrivial dynamics which become transparent on the level of the disorder quantum master equation.
We parametrize the random Hamiltonian ensemble by with a single, dimensionless disorder parameter λ, and the corresponding probability distribution p(λ). Note that, for simplicity and without loss of generality, the single parameter λ here describes, in contrast to (23), the disorder in terms of the energy level difference. In the eigenbasis {|1 , |2 } of the Pauli matrix σ z and with the initial conditions ρ jk (0) ≡ j|ρ(0)|k , j, k = 1, 2, the ensemble average dynamics are given by where φ * denotes the complex conjugate of φ. The off-diagonal dynamics are captured by the characteristic function [36] φ(t) of the probability distribution p(λ), Using (24), we can immediately derive the dynamical matrix F defined in Eq. (13). With the indices ordered as (11), (12), (21), (22), such that Time derivative and inverse of F (t) are thus easily computed and the disorder master equation (14) is where we must restrict to times t ≥ 0 (because the initial state is at t = 0). In order to obtain the desired Lindblad form, we use the Hermitian operator basis (17) with A 0 = 1/ √ 2 and A j = σ j / √ 2, with j = 1, 2, 3 and σ j the Pauli matrices. We then obtain for the matrix C defined in (19) The final form of the master equation follows directly from Eq. (21), and noting that C mn (m, n > 0) is already diagonal. One obtainṡ with the energy function and the in general time-dependent decoherence rate with ln(·) the principal branch of the complex logarithm. We thus obtain, for an arbitrary probability distribution p(λ), a dephasing master equation (28), where the energy function ϕ(t) and the decoherence rate γ(t) depend on the specific choice of p(λ). By dephasing we understand that the diagonal elements of the density matrix are timeindependent, while the off-diagonal elements evolve nonunitarily; here according to In the simplest case, the off-diagonal elements will decay monotonously, but, as we will see in the examples below, revivals can also occur. The master equation (28) thus comprehensively captures a dynamical effect which is familiar, for instance, from nuclear magnetic resonance experiments with spin 1/2 nuclei, where spatial inhomogeneities in the external magnetic field give rise to a distribution of precession frequencies. The ensemble average then amounts to averaging over these frequencies, resulting in the described dephasing given the frequencies are not all commensurate to one another. This dephasing is there characterized by the decay time T 2 [34].
Alternatively to Eqs. (29) and (30), one can express the energy function ϕ(t) and the decoherence rate γ(t) in terms of the cumulants κ (n) of the characteristic function φ(ω 0 t) (see also Appendix A): The energy function ϕ(t) depends only on the odd cumulants, which encode the degree of (a-)symmetry in the distribution (e.g. κ (1) equals the average value, and κ (3) is proportional to the skewness, i.e., "degree of asymmetry", of the distribution). The decoherence rate γ(t) is a function of only the even cumulants, which characterize the broadness of the distribution (e.g. κ (2) and κ (4) are proportional to the variance and to the kurtosis, i.e., "degree of peakedness", of the distribution, respectively). We thus learn that the coherent part of the dynamics is governed by the symmetry properties of the eigenvalue distribution, while the disorder broadness (or its strength) controls the incoherent dephasing dynamics. For example, if the disorder stems from uncontrolled parameter variations of experimental apparatuses over repeated measurements, its asymmetries would give rise to systematic deviations from the desired (coherent) dynamics, while the incurred statistical error results in the decoherence. For distributions p(λ) that are symmetric with respect to their average value, such as a Gaussian or a uniform distribution, odd cumulants of order n > 1 vanish. In that case, the coherent part of the master equation (28) is driven by the average Hamiltonian H. Irrespective of that, there will also be an incoherent contribution from the even cumulants. Therefore, the ensemble average dynamics does, also for symmetric distributions, not coincide with the unitary dynamics induced by the average Hamiltonian, i.e., exp(−iH λ t/ )ρ 0 exp(iH λ t/ ) = exp(−iHt/ )ρ 0 exp(iHt/ ).
Before evaluating the disorder-induced dephasing master equation (28) for various paradigmatic disorder distributions p(λ), let us consider its short-time approximation. The leading-order expansions of the energy function ϕ(t), Eq. (32), and the decoherence rate γ(t), Eq. (33), for short times t 1/ω 0 , yield where we used that κ (1) = λ and κ (2) This can be compared to the short-time master equation (9). With the average Hamiltonian H = λ (ω 0 /2)σ z , the short-time Lindblad operators and rates evaluate according to Eq. (10) as One can now explicitly perform the disorder integral in Eq. (9), yielding the short-time master equation for qubit spectral disorder, As expected, this coincides with the leading-order expansion (34) of the exact master equation (28), showing that the short-time master equation corresponds to the leading order of the cumulant expansion of the energy function and of the decoherence rate.
In the following, we illustrate the variety of the singlequbit dephasing dynamics that can arise from spectral disorder with four paradigmatic disorder distributions p(λ): a Cauchy-Lorentz, a Gaussian, a uniform box, and a Lévy distribution. We will see that the distinctive properties of these distributions are clearly reflected in the temporal evolution of the associated respective decoherence rates. The examples are chosen to illustrate the range of possible relations between the properties of the disorder distribution and the master equation describing the resulting ensemble average dynamics. Apart from their paradigmatic relevance, these distributions also describe realistic physical situations. For instance, as argued above, spectral disorder can account for part of the inhomogeneous broadening of the lineshapes in spectroscopy experiments. In such a setting, the Cauchy-Lorentz distribution would be related to static field inhomogeneities [34], and the Gaussian distribution to Doppler broadening or to the broadening due to impurities [37,38]. A uniform box distribution is used to model a finite-bandwidth of the distribution. The Lévy distribution is characteristic for broad-band disorder and has been used to study the spectral lines of excitonic states in molecular systems or assemblies of Rydberg atoms [39,40].

A. Cauchy-Lorentz distribution
We begin with a Cauchy-Lorentz distribution, which is familiar from resonance phenomena in the frequency domain and results in constant decay rates in the time domain, leading to an exponential decay behavior. As we will show, this expectation is confirmed in our context as well. We parametrize the Cauchy-Lorentz distribution p CL (λ) by with the (dimensionless) scale parameter σ and the (dimensionless) location parameter λ 0 . According to (25), this gives rise to the characteristic function φ CL ( The energy function (29) and the decoherence rate (30) read As we can see, the location parameter λ 0 specifies the time-independent energy function and, thus, fixes the Larmor frequency ω 0 λ 0 of the coherent precession about the z-axis induced by the Hamiltonian part of the master equation Eq.(28) (see also Fig. 5 of Appendix B). The disorder scale parameter σ sets the time-independent dephasing rate γ CL and, thereby, determines the strength of the dephasing process. More precisely, the timeindependent rate γ CL leads to a purely exponential decay with a rate proportional to σ of the off-diagonal terms of the density matrix (see Figure 2), since solving the disorder master equation (28) (31)). In nuclear magnetic resonance spectroscopy, 1/γ CL would characterize the contribution to the dephasing time T * 2 coming from the magnetic field inhomogeneity. The Cauchy-Lorentz distribution is distinguished in that it is the only disorder distribution which gives rise to both a time-independent decoherence rate γ CL and a timeindependent energy function ϕ CL , resulting in a purely Markovian Lindblad master equation obeying the semigroup property [21,22].
We note that all moments and cumulants κ (n) of the Cauchy-Lorentz distribution p CL (λ) are diverging. Nevertheless, the characteristic function φ CL (ω 0 t) of p CL (λ) (i.e. its Fourier transform) is well defined, as is it the case for any probability distribution, and consequently so are also the decoherence rate and the energy function, which are computed analytically. The energy function ϕ CL is driven by the central Larmor frequency ω 0 λ 0 because of the symmetry of the Cauchy-Lorenty distribution.
As a consequence of the non-existing moments and cumulants of the Cauchy-Lorentz distribution, the emerging dynamics exhibit a short-time behavior which differs from the one described by the short-time master equation (36), in that (36) predicts wrongly (at short times) a linearly-in-time increasing decoherence rate instead of the correct time-independent γ CL (such linear increase would induce a Gaussian rather than an exponential decay of the coherences). The correspondence to the shorttime description (36) can be restored by introducing a frequency cut-off that suppresses the algebraically decaying tails of p CL (λ), e.g., λ ∈ [−a, a], a > 0. This then corresponds to a more physical picture than the pure exponential decay, as can for example be seen in the initial Gaussian shape of the free induction decay signal in nuclear magnetic resonance experiments [34].

B. Gaussian distribution
Next, we consider the ubiquitous Gaussian (or normal) distribution p G (λ), which, in contrast to many other contexts, gives here rise to intriguing consequences. It is defined as with the (dimensionless) mean λ 0 and the (dimensionless) width σ. The energy function (29) and the decoherence rate (30) are computed from the characteristic function φ G (ω 0 t) = exp[iω 0 λ 0 − 1 2 (σω 0 t) 2 ]: The energy function ϕ G is, as for the Cauchy-Lorentz distribution, time-independent and fixed by the average value λ 0 , which gives rise to a constant Larmor precession frequency ω 0 λ 0 (see Fig. 5 of Appendix B). This is a direct consequence of the symmetry of p G (λ) with respect to its mean. In addition, the Gaussian shape of the distribution brings about a positive, linearly in time increasing decoherence rate γ G (t) whose slope scales quadratically with the variance σ. We remark that such an above all bounds increasing decoherence rate would be considered unnatural from an open-system perspective, although it can be realized with suitable engineering of the environment [41]; here, however, it arises as a natural and unavoidable consequence of the Gaussian ensemble average. On the level of the resulting average dynamics, the decoherence rate γ G leads to a fast Gaussian decay of the coherences scaling with the variance squared, 1|ρ(t)|2 ∝ exp −1/2(ω 0 σ) 2 t 2 (Fig. 2). Moreover, the positivity γ G (t) ≥ 0 of the rate for all times t ≥ 0 implies a strictly monotonous decay of the purity Tr ρ 2 (t) [20]. Furthermore, for the Gaussian distribution all cumulants of order n > 2 vanish, κ (n) = 0. Hence, the shorttime disorder master equation (36) is exact for all times, because it captures the full dependence of the ensemble average dynamics upon the first two cumulants. From the Marcinkiewicz theorem [42] it follows that the Gaussian distribution is the only possible distribution which has a finite number of non-zero cumulants, and thus it is the only distribution for which the short-time master equation (36) remains exact at all times.

C. Uniform box distribution
The uniform box distribution is illustrative for continuous distributions with a finite support (cut-off). As we will see, such distributions with a cut-off naturally induce drastic dynamical consequences such as revivals, e.g., of coherences, as it becomes manifest by the timedependence of the decoherence rates. The uniform box distribution p B (λ) is given by with the (dimensionless) location paramer λ 0 and the (dimensionless) scale parameter σ.
The corresponding characteristic function is φ B (ω 0 t) = sinc((σω 0 t)/2) exp[iλ 0 ω 0 t]. For the energy function (29) and the decoherence rate (30) one thus obtains As we can see, the energy function ϕ B is once more timeindependent and determined by the mean value λ 0 , which is due to the fact that the uniform box distribution is symmetric with respect to it. The coherent part of the dynamics is thus again a precession about the z-axis with constant Larmor frequency ω 0 λ 0 . However, in contrast to the above examples of Cauchy-Lorentz and Gaussian distributions, the width σ now not only defines a dephasing rate, but also specifies the times of the singularities τ n = 2nπ/(σω 0 ), n ∈ N + , of the decoherence rate γ B (t) arising from the box distribution (see Fig.2), the consequences of which shall be discussed below.
While the temporal behavior of the decoherence rate γ B (t) may be surprising, we emphasize that it is a direct consequence of the specific form of p B and it particularly clearly reflects the underlying ensemble average dynamics. As is evident from Figure 2, whenever the decoherence rate diverges and jumps from positive to negative values, there is a transition from decaying to increasing coherences, i.e., revivals (see also Fig. 5 of Appendix B). At the times of divergence τ n , the inverse of the dynamical map Λ −1 t is not defined. This is because, at these times, the off-diagonal terms of the ensemble average state vanish exactly. Therefore, since the master equation is a first order differential equation, a diverging decoherence rate is required to induce a revival. The divergences are formally unproblematic, because the set of diverging points τ n has measure zero and exp[−2 t 0 γ B (t)] = sinc((σω 0 t)/2) remains finite for all times t ≥ 0. Thus, the solutions to the master equation remain bounded.
The coherence revivals described by γ B (t) follow from the finite energy distribution width ω 0 σ; it has been known for long time that finite energy scales in quantum systems lead to recurrences [43,44]. Thus, any disorder distribution with a finite or closed support gives rise to revivals. On top of the revivals, the finite energy width σ ω 0 leads to a slow overall 1|ρ(t)|2 ∝ t −1 decay of the coherences.
Note that, despite the fact that γ B (t) < 0 in some time intervals, the dynamics are physical at all times, by construction from the ensemble average. This is highlighted by the fact that t 0 γ B (t )dt > 0 ∀t > 0, which can be shown to guarantee the complete positivity of the resulting dynamics using, e.g., the Jamilkowsky isomorphism [32]. The periodic divergences of the decoherence rate indicate that the dynamical map is not completely divisible. Both the negativity of the decoherence rate and the non-divisibility of the map reflect the occurrence of coherence revivals, a physical hallmark of non-Markovianity [24,45].
We remark that the short-time master equation (36) captures fully the coherent part of the average dynamics, as the exact energy function ϕ B is time-independent. However, the decoherence rate γ B (t) is approximated to first-order in time (cf. Eq. (34)). By direct comparison of the first and third order (the second order vanishes) of the cumulant expansion (33) of γ(t) in time, the range of validity of the short-time master equation is found to be restricted to times t |2γ B (t)/ ... γ B (t)| 1/2 ∝ 1/(ω 0 σ). In other words, the validity range of the short-time master equation scales inversely proportional to the disorder strength. In particular, the short-time master-equation does not capture the revivals of the coherences at times τ n ∝ n/(ω 0 σ), nor their subsequent algebraic decay.

D. Lévy distribution
The Lévy distribution stands for distributions which may be asymmetric with respect to their median value, and which have slowly decaying tails. Such distributions have been proposed to describe the effective disorder mediated by slow changes in the structured background in molecular systems [40] and are generally abundant in complex systems [46]. The Lévy distribution p Le (λ) and its characteristic function φ Le (ω 0 t) are parametrized as ϕ Le (t) and the decoherence rate γ Le (t) then follow as The location parameter λ 0 once more defines the timeindependent term of the coherent part of the master equation ϕ Le . On the other hand, the scale parameter, or more precisely its square-root √ σ, here not only fixes the rate of growth of the incoherent part γ Le (t), but also sets a rate of growth of the time-dependent term of the coherent part ϕ Le (t). This is, because, in contrast to the Cauchy-Lorentz, the Gaussian, and the uniform box distribution, the Lévy distribution is asymmetric, thus having non-vanishing odd cumulants. As a consequence, the energy function ϕ Le (t) is time-dependent. The latter leads to an also time-dependent Larmor precession frequency (see Appendix B). Figuratively speaking, thinking of spins 1/2 in a magnetic field, the asymmetry of the Lévy distribution (43) implies a statistical overweight of fast over slowly rotating spins, which induces on average a coherent rotation with a time-dependent frequency. The latter is mediated by the time-dependence of the energy function. On top of the coherent dynamics, the Lévy distribution decoherence rate leads to an overall decay 1|ρ(t)|2 ∝ exp − √ t , which is at small times faster than the Gaussian decay resulting from the Gaussian distribution, and at large times slower than the exponential decay resulting from the Cauchy-Lorentz distribution (see Fig. 2). Since γ Le (t) > 0 for all times t, the purity of the state Tr ρ 2 (t) decays strictly monotonously.
Similarly to the Cauchy-Lorentz distribution, the Lévy distribution does not possess moments and cumulants. This explains why the short-time master equation (36), predicting a time-independent energy function and a linearly-in-time increasing decoherence rate, does not capture correctly the 1/ √ t time-dependence of ϕ Le (t) and γ Le (t).

VI. SPECTRAL DISORDER IN d DIMENSIONS
We now consider the general, d-dimensional case of spectral disorder (with identical eigenvectors for all realizations, see Sec. V above). Therewith, it is for example possible to describe higher-dimensional spin states, or composite systems such as an ensemble of two-level atoms [47]. As we shall see, the particular structure of the spectral disorder directly translates into the form of the Lindblad operators.
The The disorder is then fully characterized by the eigenvalue distribution p λ = p(λ 1 , ..., λ d ) ( λ ≡ (λ 1 , ..., λ d )). Analogously to the single-qubit case, in the common eigenbasis {|j } d j=1 the diagonal terms of the ensemble average density matrix ρ(t) are time-independent, whereas the off-diagonal terms evolve according to j|ρ(t)|k = j exp[−i(λ j − λ k )ω 0 t] k , k = j. Hence, the average dynamical matrix F (t), Eq. (13), is given by with the characteristic function of the level-spacing dis- The representational matrix Q(t) of the master equation, Eq. (14), consequently reads and the non-Lindblad form of the disorder master equation, Eq. (15), is given bẏ with Π j ≡ |j j| the projectors onto the common eigenvectors of the Hamiltonians H λ . Before deriving the final Lindblad form, it is convenient to compare the short-time expansion of (49) with the short-time disorder master equation (9). The leadingorder in time approximation of (d/dt) ln(φ * jk (t)) in (49) reads Inserting (50) into (49), one can deduce the non-diagonal Lindblad form for the short-time approximation: On the other hand, with the ensemble-averaged Hamiltonian H = ω 0 d j=1 λ j |j j|, one obtains for the shorttime Lindblad operators and decay rates (10) Inserting (52) into (9) and performing the integrals over the disordered eigenvalues, we find that, as expected, the short-time master equation coincides with the leadingorder expansion (51) of the master equation.
We proceed now towards the general Lindblad form of the spectral disorder master equation. Following the matrix approach outlined in Section IV, we obtain in terms of a Hermitian, traceless operator basis {A j } d 2 −1 j=0 the effective Hamiltonian H(t) = (i /2) C (t) −C † (t) (cf. Eq. (21)) with Since φ jk = φ * kj and A m,jj ∈ R, the decoherence matrix Γ(t) = Γ † (t) is, as expected, Hermitian and can therefore be diagonalized for given characteristic functions φ jk (ω 0 t). The resulting Lindblad operators in the diagonal form of the master equation are then linear combinations of the G l d and thus also diagonal operators. Consequently, general spectral disorder always results in dephasing dynamics in the state basis {|j } d j=1 defined by the disorder ensemble. It is thus a natural generalization of the single qubit with spectral disorder considered in the previous section. The populations are time-independent, whereas the off-diagonal terms of the density matrix evolve non-unitarily. This result can be understood as follows: For each single realization of the disorder, the populations of the individual Hamiltonians' eigenstates are time-independent and thus also remain time-independent after the averaging, whereas the oscillating off-diagonal terms undergo dephasing due to the averaging over the disorder in the eigenvalues.
Although there is no direct cross-talk between the coherences on the level of the master equation, their dynamics may exhibit strong correlations, mediated by the correlations among the elements of the decoherence matrix. Furthermore, since the decoherence matrix depends only on the level-spacing distribution, all disorder distributions p λ which give rise to the same level spacing statisticsp jk (∆λ) induce the same disorder master equation.
In the following, we have a more detailed look at two specific types of spectral disorder, namely global spectral disorder, which is characterized by a single disorder parameter, and fully uncorrelated spectral disorder.

A. Global spectral disorder
Global spectral disorder represents a direct generalization of the previously studied spectrally disordered single qubit to d dimensions, in the sense that the spectrum is only subject to a single disorder parameter. This may for example describe an array of N uncoupled two-level atoms in a common, static magnetic field, which varies slowly in intensity over time, giving rise to collective dephasing [47].
We consider the ensemble of random Hamiltonians with {|j } d j=1 the common eigenbasis and ω 0 j the eigenvalues of the reference Hamiltonian H 0 = d j=1 ω 0 j |j j|. The single disorder parameter λ is distributed according to p(λ). Note that, for convenience, the global disorder definition (56) differs from the definition (45) for general spectral disorder. The two definitions correspond via the identification ω 0 λ j = ω 0 j λ, which relates the probability distributions via p λ = {p(λ) if λ j = λ(ω 0 j /ω 0 ); 0 otherwise}. In terms of the more natural characterization (56) of global spectral disorder, the dynamics are characterized by a single characteristic function φ evaluated at the reference Larmor frequency differences ω 0 jk ≡ ω 0 j − ω 0 k , which now mediate the dependence upon the indices j, k, All previous results for general spectral disorder can now be taken over by identifying the characteristic functions φ jk (ω 0 t), (47), with φ(ω 0 jk t), (57). As a consequence of the single joint disorder parameter for all eigenvalues, the dynamics of the coherences are strongly correlated in terms of the characteristic function (57), j|ρ(t)|k ∝ φ * (ω 0 jk t). On the other hand, they also depend on the Larmor frequency difference ω 0 jk connecting the levels j and k, which renormalizes their decoherence time. Therefore, all coherences will decay with the same envelop function, but at different velocities.
Interestingly, if there are degeneracies in the spectrum of the reference Hamiltonian H 0 , i.e., at least two eigenvalues ω 0 j and ω 0 k coincide, the off-diagonal terms of the density matrix ρ(t) associated with these degenerate eigenvalues are completely time-independent (φ(ω 0 jk t) = 0) and, hence, survive the averaging procedure, even in the asymptotic limit. This mechanism can for example be exploited in order to find long-lived entangled states [47].
Finally, it is obivous that the characteristic function (57) is the same as the characteristic function (25) for a single qubit with spectral disorder, given the probability distributions coincide. Thereby, the Hamiltonian H(t) and the decoherence matrix Γ mn (t) are for global spectral disorder directly obtained by identifying the real and imaginary part of − d dt ln[φ * (ω 0 jk t)] in Eq. (53) and Eq. (54) with the qubit decoherence rate 2γ jk (t), Eq. (30) and the qubit energy function (2/ )ϕ jk (t), Eq. (29), respectively, where the indices j, k indicate that ω 0 is replaced by ω 0 jk . In particular, for a Gaussian distribution p G (λ) with width σ, cf. Eq. (39), the qubit rate and energy function (40) One can then easily derive a diagonal Lindblad master equatioṅ where the single Lindbald operator is given by the reference Hamiltonian H 0 . Since the global spectral disorder master equation for a Gaussian distribution, Eq. (58), is at most linear in time, it coincides with the corresponding short-time master equation obtained from (51). We thus find that in this case the short-time dynamics remains valid at arbitrary times. Interestingly, due to the absence of higher cumulants, the master equation can be generalized to time-dependent noise λ(t), as was studied in [48].

B. Uncorrelated spectral disorder
We finally consider a spectral disorder ensemble (45) with a fully uncorrelated eigenvalue distribution, where the probability distributions p j , j = 1, . . . , d, may in general differ from one another. The characteristic function (47) then factorizes and The characteristic function φ j (ω 0 t) specifies the Hamiltonian via (53) and the decoherence matrix via (54). We can then derive a diagonal Lindblad form of the master equation, where Im d dt ln φ * j (ω 0 t) Π j and the Π j = |j j| are again the projectors onto the common eigenvectors of the Hamiltonians in the ensemble. Since there are no correlations in the considered spectral disorder, each decoherence rate is equal to γ j (t) ≡ −2Re d dt ln φ * j (ω 0 t) and thus depends only on the characteristic function corresponding to the j th eigenvalue. Hence, as opposed to the previous strongly correlated case of global spectral disorder, the different coherences will in general exhibit disparate decay patterns. Note that the above diagonal Lindblad form (60) can also be directly derived from the projector form (49).

VII. UNITARILY INVARIANT DISORDER
In the previous section, we studied the ensemble average dynamics resulting from spectral disorder, where all Hamiltonians in the ensemble shared a common eigenbasis. We now generalize to the case that also the eigenstates are subject to disorder. For example, think again of the ensemble of spins 1/2 in a static magnetic field from Section V, but now in addition to the magnitude also the orientation of the field is random. The latter can for instance be induced by local impurities. More generally speaking, disorder in the eigenstates is implied whenever the Hamiltonians in the ensemble (1) do not all commute with one another. This is, for instance, typical for systems with a kinetic hopping term and a disordered on-site potential, such as, e.g., the Anderson model [49]. In the following we elaborate how disorder, not only in the eigenvalues, but also in the eigenvectors, can affect the structure of the master equation describing the ensemble average. Concretely, we focus on an isotropically randomized eigenstate distribution.
Let us consider an ensemble of d×d Hamiltonians with uniformly distributed eigenvectors and with eigenvalues distributed independently of the eigenvectors, where D λ = ω 0 diag(λ 1 , ..., λ d ) and W W † = 1. The isotropic eigenvector distribution is reflected by the invariance of the probability distribution p λ with respect to the unitary transformations W , p 1, λ = p W, λ = p λ . Prominent members of the disorder type (61) are the unitarily invariant random matrix ensembles [50][51][52]. Generally speaking, random matrices have a wide range of application in the study of the statistical properties of complex quantum systems [53][54][55][56].
In order to understand the impact of the uniform eigenvector distribution on the ensemble average dynamics, we now derive the master equation for the unitarily invariant disorder ensemble. The evolution of a single realization is given by ρ W, λ (t) = W exp −(i/ )D λ t W † ρ(0)W exp (i/ )D λ t W † . The ensemble average state is then obtained by integrating both over the unitaries W in terms of the Haar measure d µ (W ), and over the eigenvalue distribution p λ = p(λ 1 , ..., λ d ): The Haar measure integral in Eq. (62) can be conducted using the two following results which are obtained using the Weingarten calculus for unitary groups [57,58], and Using these, one obtains for the time evolution of the ensemble-average state the dynamical map It describes a depolarization channel, i.e., the mixing of the initial state ρ(0) with the maximally mixed state 1/d, with the time-dependent mixing probability a(t) (a(0) = 1); the latter is determined by the eigenvalue distribution p λ in terms of the sum over the characteristic functions for all level spacings (λ j − λ k ), Eq. (47), which we denote as Note that here the bar indicates averaging with respect to the eigenvalue distribution p λ only. Following the general procedure outlined in Section IV, one obtains the non-Lindblad master equation for unitarily invariant disorder, The non-Lindblad master equation (67) is described by the representation matrix Q jkrs = (ȧ(t)/a(t))(δ jr δ ks − (1/d)δ jk δ rs ), cf. Eq. (15). In Section IV, we provided a systematic method to derive the general diagonal Lindblad form based on the representation matrix Q. In the present case, however, it is possible to deduce the diagonal Lindblad forṁ directly from Eq. (67) by requiring that the Lindblad operators satisfy j L j (t)ρ(t)L † j = 1/d and j L † j L j = 1 (i.e. the Lindblad operators form a basis in operator space). The single depolarization rate γ(t) is then given by A possible choice for the Lindblad operators L j are the Hermitian Gell-Mann matrices G j d 2 −1 j=1 [30] (see also Eq.(17)), complemented by the identity matrix and with a suitable normalization We thus find that the master equation (68) for unitarily invariant disorder does, first of all, not have a coherent Hamiltonian term. This is because averaging over a uniform distribution of eigenvectors cancels out all the terms arising from the (a-)symmetry in the eigenvalue distribution, i.e., the odd moments and cumulants, that generate the coherent dynamics (see Section V). Moreover, the Lindblad operators L j are independent of the eigenvalue distribution; in other words, they are completely specified by the uniform eigenvector distribution alone, describing the isotropic dynamics towards the maximally mixed state which is characteristic of depolarization. The single depolarization rate γ(t), on the other hand, is characterized by the eigenvalue distribution p λ , more specifically by the level-spacing distribution. Therefore, different eigenvalue spacing statistics will lead to different mixing probabilities, as we will show in the examples below.
Before, it is interesting to see how the structure of the master equation (68) is recovered in the short-time limit from the general expression (9). In the presently considered case of a unitarily invariant ensemble (61), one obtains, using (63), that the average Hamiltonian is H = d λ d µ (W )H W, λ = Tr D λ 1/d. According to (10) the short-time decay rate and the Lindblad operators are thus given by With (71), and using d µ (W )W X 1 W † X 2 = (Tr[X 1 ]/d)X 2 and (64), one can perform the Haar measure integral in Eq. (9). Based on a suitable basis change, we then obtain the short-time diagonal Lindblad master equatioṅ This short-time master equation, which is obtained with the method outlined in Section III, is in agreement with the leading-order expansion in time of the master equation (68) for t 1/ω 0 . Note that the Lindblad term of the short-time master equation (72) is composed of the same Lindblad operators L j , which are determined only by the uniform distribution of the eigenvectors, as the finite-time master equation (68). On the other hand, the decoherence rate of the short-time description in (72) contains only the first-order in time term of the exact rate (69). Therefore, differences between both equations must stem from the eigenvalue distribution specifying γ(t). Possible deviations of the short-time approximation from the exact dynamics at larger times can thus be traced back to p λ . By direct comparison of the first-order and third-order (second order vanishes) terms of the short-times expansion of γ(t), we derive that the validity of the short-time master equation is limited to times which, as we will see in the example below, may be shorter than 1/ω 0 .
In summary, for unitarily invariant disorder, the distribution of the eigenvectors fully characterizes the Lindblad operators, i.e., the character of the decoherence process, whereas the distribution of the eigenvalues determines the temporal evolution of the corresponding decoherence rates. We will now treat in more detail the examples of the Gaussian unitary ensemble (GUE), which is often employed to uncover generic features of quantum chaotic systems [59], and, for comparison, the "Poissonian ensemble" (PE), which is meant to model an integrable counterpart [50]. These ensembles are distinguished in their distribution of the eigenvalues p λ and, as we shall see, the corresponding master equations are accordingly distinguished by the time-dependence of their decoherence rates.

A. Poissonian ensemble
As a first example, we consider an uncorrelated eigenvalue distribution with each eigenvalue ω 0 λ j following a uniform box distribution p B (λ j ) of width σ (cf. Eq. (41)), This gives rise to the Poisson distributed eigenlevel spacing of the PE. To determine the depolarization rate (69), we note that for permutation symmetric distributions such as (74) one has where is the two-point correlation function of the distribution p λ [52]. For the Poissonian ensemble with the probability distribution (74), the two-point correlation function can be evaluated analytically and reads R 2 (λ 1 , λ 2 ) = d(d − 1)p B (λ 1 )p B (λ 2 ). One then obtains the depolarization rate (69) In the limit of large dimensions, d 1, the rate converges to which may be compared to the single-qubit decay rate (42) for uniformly distributed spectral disorder. In Fig. 3, we show how for different finite dimensions d the depolarization rate decays algebraically while oscillating around the origin with a period T = (2π)/(ω 0 σ). In the highdimensional limit, we reobtain the periodic divergences already encountered in the case of a single qubit subject to uniform-box-distributed spectral disorder. This underlines that exactly at the times when the rate diverges all coherences completely vanish before undergoing a subsequent revival (see Fig. 3). Note that the general validity condition (73) of the short-time master equation evaluates for the PE to t 1/(ω 0 σ), scaling inversely proportional to the disorder strength.
The dynamical properties of the ensemble average state for unitarily invariant disorder are captured best in terms of its purity, We remind the reader that, by definition, the purity lies in the interval [1/d, 1], where the lower bound corresponds to the maximally mixed state, and the upper to a pure state. In case of the Poissonian ensemble, the mixing probability, Eq. (66), is a P (t) = (1 + d sinc 2 ω0σ 2 t )/(1 + d). As we can see in the insets of Fig. 3, each time the depolarization rate turns negative a revival of purity occurs. In addition, one finds that for finite dimensions the asymptotic state is not the maximally mixed state (Tr ρ 2 (t) = 1/d), in spite of the tendency of the depolarization dynamics to drive any initial state towards the maximally mixed state. Formally speaking, this is due to the degenerate terms in χ(t) which lead to the constant factor 1/d in (75). Consequently, the larger the dimension, and hence the more "depolarization directions" there are, the smaller the contribution 1/d in (75) gets, and the closer the state is driven towards the maximally mixed state. Being more precise, the asymptotic limit of the purity is given by  Tr ρ(t) 2 |, showing that the purity oscillations are directly correlated with the oscillations of the rate γP(t) around zero. We chose the uniform box distribution width σ = 4 for comparison with the GUE as in [60].

B. Gaussian unitary ensemble
We consider as a second example the Gaussian unitary random matrix ensemble. The joint eigenvalue probability distribution for the GUE [51] reads where A is a normalization constant fixing the average density of states at small energies. Following the conventions in [57,60], we set the normalization so that H 2 jk /( ω 0 ) 2 = 1 d (thereby fixing the variance). Since the Gaussian distribution (81) is invariant with respect to permutations of the eigenvalues λ j , we follow the same procedure as in the previous section in order to compute the sum χ(t) of the level-spacing characteristic functions, cf. (75). The two-point correlation function for the GUE can be expressed as [57,60] where ϕ n (x) = (2 n n! 2π/d) −1/2 e −x 2 /2 H n ( d 2 x) are the harmonic oscillator eigenfunctions, with the Hermite polynomials H n (x) ≡ exp x 2 (−d/dx) n exp −x 2 . Inserting R 2 (λ 1 , λ 2 ) into Eq. (75), one can evaluate χ(t) for finite dimensions d. This then allows one to compute the depolarization rate γ GUE (t), Eq. (69), and the GUE mixing probability a GUE (t), Eq. (66), specifying the dynamics of the state purity (79). For example, for d = 2 the depolarization rate reads For the high-dimensional limit, we use the approxi- denotes the Bessel function of the first kind, in order to evaluate the GUE depolarization rate γ GUE (t) and the mixing probability a GUE (t).
In Fig. 4, we show the time evolution of the GUE depolarization rate γ GUE (t) for dimensions d = 2, 4, 8, ∞. For finite dimensions, γ GUE (t) first exhibits algebraically damped oscillations, before approaching the x-axis exponentially from below. In the limit of infinite dimensions, the GUE rate follows a similar temporal evolution as the PE rate, Eq. (78), including periodic divergences. This reflects a universal property of random matrix ensembles: in the limit of high dimensions, the eigenvalue statistics only depend on the density of states and the type of symmetries fulfilled by the random matrices [61] (in the present case, the unitary invariance).
In Fig. 4 we also show the time evolution of the purity (79). For finite dimensions, the purity initially decays to the minimum 1/d of the maximally mixed state, in contrast to the Poissonian case (see Fig. 3). Afterwards, it increases again towards its asymptotic value (80). The generic purity revivals at short times, which are clearly visible in the insets of Fig. 4, are directly related to the oscillations of the depolarization rate γ GUE (t) (for more details see Appendix C). Additionally, we find that the purities approach their asymptotic values exponentially, free of oscillations, from the point when the corresponding decoherence rates decay exponentially to zero. In the high-dimensional limit, we recover a similar behavior as for the Poissonian ensemble, including an infinite number of purity revivals and an algebraic decay towards the asymptotic state. Finally, the asymptotic purities of the GUE coincide, because of the absolute integrability of the eigenvalue distribution of the GUE (81) and due to the Riemann-Lebesgue Lemma, for all dimensions with the values for the PE (c.f. Eq. (80)).
In summary, systems subject to unitarily invariant disorder all lead to the same asymptotic purity, but follow-  For finite dimensions, the rate undergoes a finite number of oscillations and a subsequent exponential decay to zero from below. Bottom panel: Decay of the purity Tr ρ 2 (t) for a pure initial state, Tr ρ 2 (0) = 1. For finite dimensions, the asymptotic values of the purities correspond to Eq. (80), which is larger than 1/d (indicated by horizontal, black, dotted lines). Inset: Double-logarithmic plot of the approach to the asymptotic state Tr ρ 2 (t) − lim t→∞ Tr ρ 2 (t) . The purity oscillations are directly correlated with the oscillations of the rates γGUE(t) (dotted, vertical, gray lines). Note that, for finite dimensions, the peaks are shifted w.r.t. their corresponding rates γGUE(t) in the top inset, because of the subtraction of the asymptotic purity value.
ing different transient paths characterized by the different spectral disorder distributions. For finite dimensions, the regular Poissonian ensemble leads to an overall algebraic decay towards the maximally mixed state, in contrast to the exponential decay of the chaotic Gaussian ensemble. Similar behavior was noted in the semiclassical description of the decay of autocorrelation functions of chaotic and regular quantum systems [62].

VIII. CONCLUSIONS
We introduced quantum master equations as natural and viable means to conceptually grasp the incoherent effective dynamics of quantum systems described by an ensemble of time-independent, disordered Hamiltonians. To this end, we presented a general method to derive these master equations as they emerge from the disorder average. Therewith, we provided a new tool to study the average dynamics of disordered quantum systems both on transient and asymptotic time scales. Our approach furthermore naturally separates coherent and incoherent components of the dynamics, and associates characteristic time and coupling constants of the effective time evolution with the statistical features of the generating disorder.
In the paradigmatic case of spectral disorder, we found that the common basis of eigenstates of the ensemble Hamiltonians specifies the decoherence process, i.e., dephasing, while the eigenvalue distributions determine the corresponding time-dependent decoherence rates. Similarly, in the illustrative case of unitarily invariant disorder, we obtained that the isotropic eigenstate distribution leads to a depolarization process with a rate again given by the eigenvalue distribution. These two cases were solved exactly. It is clear that treating more complex and realistic disorder models will require some degree of approximation. For example, the master-equation in Lindblad form presented in Section III, which is valid in the limit of short times, has been used in [17] to study homogeneous disorder models. In order to extend its validity beyond short times, one could employ a Born-Markov-type approximation, as done in [63] to derive an evolution equation for diffusive spin-transport in a disordered medium. Another perspective to describe complex disordered system may consist in complementing our master equation approach with techniques from disorder transport theory such as perturbation theory or diagrammatic methods as used, e.g., in the study of weak localization [64]. Nevertheless, the generic examples here treated in detail clearly illustrate that the master equation description of the ensemble-averaged dynamics is not only a new approach to the dynamics of disordered systems, but also provides useful and non-trivial information, e.g., on the coherent and incoherent contributions to the effective dynamics, respectively. The latter cannot be readily extracted from the asymptotic, disorder-averaged state alone.
From a more general point of view, we conclude that the character of the decoherence process is determined by the structure of the disorder, i.e., the eigenstate distribution, while the temporal course of the master equation follows from the eigenvalue distribution. The latter then mediate dynamical features such as revivals. Therefore, knowing these master equations, we are able to characterize concisely the collective dynamics of the ensemble average. Conversely, their experimental inference may allow one to identify the underlying disorder properties. In other words, the experimental monitoring of the effective dynamics on transient time scales could be used as a diagnostic tool to uncover the characteristic features of the disorder with the help of the master equations introduced in this article.
The necessity to average over an ensemble of disordered Hamiltonians constitutes a ubiquitous issue in the quantum theory of complex systems, be it for experimental or for conceptual reasons, and to this end the presented master equation description may be decisive in order to obtain a thorough understanding of the associated dynamics. Besides the discussed occurrence of spectral disorder in trapped ion systems exposed to fluctuating external magnetic fields [35], we would like to mention, as other interesting experimental examples, the transport of electrons or excitations in biomolecular systems subject to slow conformational changes [55,65], cold atom implementations of quantum walks, where residual thermal fluctuations lead to an effective differential light shift [66], and the propagation of photons in the presence of a turbulent atmosphere, giving rise to decoherence in the orbital angular momentum degree of freedom [67,68]. On the conceptual side, we mention the characterization of the ensemble average dynamics in the seminal Anderson model [49], coherent backscattering [69,70], and the identification of robust features in boson sampling systems [71] based on random matrices [72].
Besides understanding the dynamical impact of averaging over an ensemble, the master equation approach also suggests that disorder can serve as a useful resource to generate quantum dynamical maps of interest. Each disorder master equation gives rise to a family of random unitary maps, which are important tools in quantum information theory [18,26,[73][74][75]. More generally speaking, it may be insightful to clarify the relationship between the ensemble average dynamics of disordered systems and the reduced dynamics of open systems, i.e., to identify environmental models giving rise to the same system dynamics as the ensemble average. This will help to clarify whether encountered dynamical properties must be traced back to disorder or may also be a consequence of the coupling to an environment. Our approach thus establishes a non-trivial connection between the theory of disordered and open quantum systems, and of quantum information.

Appendix D: Equivalence of time-local and non-local form
We present the essential steps showing that master equations in local form (c.f. Eq. (6)) and master equations with a memory kernel are equivalent to one another [23,25].
Suppose the dynamical map ρ(t) = Λ t−t0 [ρ(t 0 )] is generated by an integro-differential master equation with memory kernel K s (t − t 0 ) Using the formal inverse of the dynamical map, Λ −1 t−t0 , we obtain with the mapΛ t−t0 = t t0 ds K s (t − t 0 ) • Λ s being a function of t − t 0 only (and not of s anymore). The last equality (D5) is precisely our time-local master equation (11). Therefore, the intergro-differential form and the time-local form are equivalent, and the latter is perfectly able to describe memory effects. This is essentially because the time-local generator explicitly depends on the time-difference t − t 0 and not only on t (this can easily be overlooked as one typically sets t 0 = 0). On a technical level, going from one form to another can for example be done by taking the Laplace transform.