Protecting quantum coherences from static noise and disorder

Quantum coherences are paramount resources for applications, such as quantum-enhanced light-harvesting or quantum computing, which are fragile against environmental noise. We here derive generalized quantum master equations using perturbation theory in order to describe the effective ensemble-averaged time-evolution of finite-size quantum systems subject to static noise on all time scales. We then analyse the time-evolution of the coherences under energy broadening noise in a variety of systems characterized by both short and long-range interactions, by strongly correlated and fully uncorrelated noise -- a single qubit, a lattice model with on-site disorder and a potential ladder, and bosons in a double-well potential with random interaction strength -- and show that couplings can partially protect the system from the ensemble-averaging induced loss of coherence. Our work suggests that suitably tuned couplings could be employed to counteract part of the dephasing detrimental to quantum applications. Conversely, tailored noise distributions can be utilized to reach target non-equilibrium quantum states.


I. INTRODUCTION
For applications such as quantum-enhanced efficient light-harvesting [1] or quantum computing [2] coherences are essential resources that in general are fragile against environmental noise [3].Among these, the static noise sources (often also termed disorder) are typically seen as technological imperfections, whose impact on the effective dynamics of the coherences can however be dramatic [4].The latter can also be actively exploited, for instance for random lasing [5], channel coding of a measurement [6], or novel multi-component materials with unprecedented properties [7].Hence, a deep understanding of the effects of the static noise on the coherences, and in particular their time-evolution, is crucial for the development of quantum technologies [3].
Characterizing the dynamics of noisy systems requires a statistical approach, as single realizations cannot give reliable predictions on reproducible features of the system as a whole.In the case of static noise, one studies an ensemble of unitary trajectories each one characterized by a random time-independent Hamiltonian.The most natural and accessible quantity is the average over this ensemble, since it is often difficult or even impossible to measure single realizations for several successive experimental runs as would be required to access higherorder correlations.As was established and extensively discussed in Ref. [4], the ensemble-averaging procedure implies an average over the accumulated phases associated with the eigenstates of every realization of the underlying random Hamiltonian which induces a loss of phase information, hence dephasing, i.e., a specific type of decoherence [2,8].In other words, the ensembleaveraged state must be described in terms of a density matrix as it evolves non-unitarily, in contrast to the state of individual realizations of the noise which evolve unitarily.A simple example of averaging-induced decoherence is the decay of the transverse-polarization in nuclear magnetic resonance experiments due to magnetic field inhomogeneities [9].Note that static noise can only give rise to unital (non-unitary) ensemble-averaged dynamics, which excludes certain types of decoherence such as amplitude damping.
An appropriate theoretical framework to describe the measurable ensemble-averaged dynamics, and in particular the evolution of the coherences (off-diagonal elements of the density matrix in a specific basis), is provided by generalized (quantum) master equations [4].We here present a master equation that describes the ensembleaverage dynamics of finite disordered quantum systems that can be treated with time-independent perturbation theory in the absence of noise.Note that the small parameter needs not be specified more precisely ; it is sufficient to assume that the system in question is well described by standard perturbation theory.Hence, our approach is quite general and allows to describe the dynamics of a large variety of noisy systems in terms of effective, physically meaningful (time-dependent) shifts and (time-dependent) rates.A key point of the approach is that the perturbative expansion is done on the level of the state's representation -more precisely, on the level of the dynamical matrix which contains the same information as the Choi matrix [10] -, and not on the level of the (Von Neumann) equation of motion for the density matrix.We thus not only obtain a description of the ensemble-averaged dynamics on transient times, but also of the asymptotic state.
We will derive perturbative master equations to study the time-evolution of the coherences in a variety of systems with static noise in the energies (diagonal noise)a broadened two-level system, a lattice model with onsite disorder, and bosons in a double-well potential with random interaction strength.We demonstrate that the coherences of the levels coupled by the perturbation are protected from complete averaging-induced decoherence (see illustration Fig. 1).As we will discuss in more details, this occurs because the ensemble-averaged dynamics describe a dephasing process in the eigenbasis of the noise-free Hamiltonian, and consequently the asymptotic state is the projection of the initial state onto the eigenbasis of this Hamiltonian, which may then exhibit nonvanishing coherences.The article is organized as follows.The static noise model and the ensemble-averaged dynamics are defined in Section II.The general perturbative master equation is derived in Section III, and the particular form for diagonal static noise in Section IV.The effects of the noise on the dynamics of the coherences is studied for various systems in Section V, and Section VI concludes.

II. STATIC NOISE AND ENSEMBLE-AVERAGED DYNAMICS
Our model considers quantum systems that can effectively be described by an ensemble of finite-size random d-dimensional Hamiltonians [11], which accounts for the descriptions of a large variety of experimental scenarios -quantum systems coupled to a classical environment [11][12][13], slow fluctuations in the experimental system parameters over time [14,15] or generic statistical distribu-tions of parameters usually referred to as disorder [16][17][18].Here we focus on the particular case that for each realisation of the noise, the system is described by a random Hamiltonian Ĥλ with an unperturbed, diagonalizable part Ĥ(0) λ , and a perturbation α Vλ , where α 1 is a dimensionless perturbation parameter and λ labels the individual realisations.The random Hamiltonian ensemble is then characterized by with p λ the probability density associated with the realization λ.The ensemble-averaged dynamics are characterized by the average density matrix where ρ(0) is the initial state (independent of the realisation of the noise), and overline • denotes ensembleaveraged quantities.Note that in this manuscript we will use the terms static noise and disorder interchangeably, as in our context they equally give rise to a statistical ensemble description of the form Eq. ( 1).
For simplicity, we further assume that the unperturbed Hamiltonians Ĥ0 λ are diagonal in the same eigenbasis {|j } d j=1 for all realisations of the noise.This assumption is for convenience, and serves to identify a specific basis in which to represent the master equation.We thus consider Hamiltonians consisting of a diagonal term Ĥ0 λ = d j=1 E 0 λ,j |j j| and a small coupling potential α Vλ = d j,k=1 V λ jk |j k|.

III. DERIVATION OF THE PERTURBATIVE MASTER EQUATION
In this section we derive a perturbative master equation characterized the ensemble-averaged dynamics, Eq. ( 2), of the random ensemble, Eq. ( 1).As a starting point we derive the (non-degenerate) perturbative expansion of the unitary dynamics [19,20] for a single realisation of the static noise and build the associated perturbative dynamical matrix [4,21].
To proceed further, we represent the unitary timeevolution superoperator U λ (t) by its dynamical matrix F λ (t) and identify the Liouville states |ρ} with vectors ρ [4], such that Concretely, the dynamical matrix F λ (t) for a single realisation of the static noise is obtained by expanding U λ (t) in an orthonormal basis and identify the matrix elements.Given {|jk}} an orthonormal basis in Liouville space, and {|j } the common eigenvectors of the unperturbed Hamiltonians Ĥ0 λ , we have The dynamical matrix has double indices (jk) with j, k = 1, . . ., d ordered as (11), (12), . . .(1d), (21), . . ., with d the dimension of the system.We can now expand the evolution operator Eq. ( 4) in the basis {|jk}} of the unperturbed Hamiltonians, and identify the corresponding dynamical matrix for each order in α: The matrix elements of each order can be evaluated separately (detailed calculations can be found in [22]).
The ensemble average of the dynamical matrix, is obtained by taking the weighted integral over all realisations of the noise of the perturbative dynamical matrix, Eq. ( 7), and reads Note that the ensemble-averaged dynamics, Eq. ( 2), are obtained from ρ (t) = F (t) • ρ(0), which in terms of density matrix elements reads j|ρ(t)|k = d r,s=1 F jk,rs (t) r|ρ(t)|s .
To derive a corresponding perturbative disorder master equation, we follow the general procedure introduced in [22], which is analogous to the time-convolutionless expansion (TCL) approach from the theory of open quantum system [8].The first step is to compute the matrix to obtain a time-convolutionless form of the generator, i.e. ˙ ρ = Q (t) • ρ .The time derivative Ḟ (t) is directly obtained by taking the time derivative of Eq. (8).Computing the inverse F −1 (t) analytically is, in general, not possible.We therefore make use of the perturbative nature of the dynamics and express the inverse with the help of the Neumann series [23] as The expression for Q (t) = Ḟ (t) • F −1 (t) then yields The matrix Q (t) fully characterizes the perturbative master equation for the ensemble-averaged dynamics, Eq. ( 2), and can be truncated at the appropriate order in the perturbation parameter α.If desired, a Lindblad form of the master equation can be obtained by expanding Q(t) in a basis of traceless, Hermitian operators [4,8,24].Alternatively, one can express the master equation in terms of the density matrix components and then collect and rearrange the terms in a suitable way.

IV. PERTURBATIVE MASTER EQUATION FOR DIAGONAL STATIC NOISE/DISORDER
In this section we derive the explicit form of the perturbative master equation to first order for d-dimensional systems with diagonal (spectral) static noise and a weaker, noise-and time-independent perturbation.More precisely, the random Hamiltonians read where α |E 0 λ,j − E 0 λ,k | for all j, k such that V jk = 0.In other words, the coupling potential between two unperturbed eigenstates is much smaller than their energy difference.The system can thus be treated with nondegenerate perturbation theory [20].
Since we assume the perturbation to be independent of the noise ( Vλ = V ), the eigenvalues E 0 λ,j of the unperturbed Hamiltonian Ĥ0 λ are the only random variables.Indeed, we remind the reader we assume the eigenvectors {|j } of Ĥ0 λ to be independent of the noise realisation (on the contrary the eigenvectors of Ĥλ are obviously modified by the noise).If the eigenvectors of Ĥ0 λ would also be subject to noise, one would have to expand the dynamical matrix, Eq. ( 7), in a noise-independent basis before taking the ensemble average.This does not change the derivation procedure but the computations would, in general, become more intricate.Furthermore, we define, without loss of generality, that in the eigenbasis of the unperturbed Hamiltonian Ĥ0 λ , the diagonal elements of the perturbation vanish, j| V |j = 0 ; possible contributions of V to the diagonal can be absorbed into the definition of the ensemble-averaged value of the unperturbed Hamiltonian Ĥ 0 = dλ p λ Ĥ0 λ .
Interestingly, and as will become clear in the next section, the first order correction to the master equation matrix Q (t) is enough to capture the main features of the ensemble-averaged dynamics at all times for random ensembles of the form of Eq. (12).The first order in α master equation can be expressed as (for a detailed derivation see App.B) with the diadic operators Πjk := |j k|.The (timedependent) rates in the equation above are functions of the phase factors and their average ϕ jk (t) = d λp λ ϕ jk (t).The rates are defined as = − Γ * kjjr (t).
Here the symmetry Γ jkrj = −Γ * kjjr is obtained by permuting the first two indices, jk → kj, and the last two indices, rj → jr, and reflects the Hermiticity of the density matrix (the dynamical matrix F (t) has the same symmetry as Γ(t)).The first-order in α part of the master equation, Eq. ( 13), is manifestly separated into a fully coherent contribution from the perturbation potential V (second term, first line), and additional incoherent terms with decoherence rates ∝ γ(t) and ∝ Γ(t) (second and third lines).
In order to better understand the dynamics arising from the master equation, Eq. ( 13), we express it in terms of the density matrix elements.For the coherences (offdiagonal elements j = k) we have The term proportional to Υ jk (t) describes the effects of pure diagonal disorder (α = 0) [4], and leads to timedependent dephasing in the eigenbasis of Ĥ0 .All other terms arise as a consequence of the perturbation potential V .The rate γ jk (t) is associated with the dynamical coupling of the coherences to the population differences, whereas the rate Γ jkrj (t) governs the second-order coupling of the coherences to the other off-diagonal terms of the density matrix.At the same time, for the populations (diagonal terms j = k) we have Hence the populations do evolve in time (as opposed to the case when α = 0), but since Eq. ( 19) does not contain any noise-dependent decoherence term, their evolution is not directly affected by the noise.They only indirectly feel the noise through their coupling to the coherences via the perturbative potential V .The ensemble-averaging for ensembles of the type (12) thus leads to complex, in general time-dependent, dephasing (decoherence) processes.

A. Statistical interpretation
The ensemble-averaged phase factors ϕ jk (t) that characterize the decoherence rates, Eqs. ( 15)-( 17), can be expressed in terms of the complex conjugate of the characteristic function [25] of the probability density distribution q jk (∆ jk ) := Hence, the time-dependent properties of the ensembleaveraged dynamics can be traced-back to the properties of the characteristic functions of the pair-wise eigenvalues differences distributions.We remark that higher order moments (involving more than two eigenvalues) are necessary if one considers contributions beyond the firstorder in the perturbation.
Conveniently, the characteristic function can be described in terms of its cumulants [26] which yield a direct physical interpretation [27] -the odd cumulants capture the (a)symmetry of the distribution, and contribute to the coherent (Hamiltonian) part of the master equation, while the even cumulants describe the width of the distribution (or strength of the noise), and characterize the decoherence rates Eqs.( 24)-( 26), i.e., the time-dependence and speed of the dephasing process.
Note that while the characteristic function describes the time-dependence of the decoherence process, the structure of the Hamiltonian, and in particular the nature of the coupling V , determines the genre of the decoherence [4].Indeed, the coupling potential determines the Lindblad operators of the master equation ( 13), i.e., the type of dephasing/decoherence dynamics.This is consistent with findings in [4] that the structure of the disorder defines the form of the Lindblad operators.

B. Range of validity
In the derivation of the general perturbative master equation (c.f.Section III) we assume the conver-gences of the Neumann series used to expand the inverse of the ensemble-average dynamical matrix F , i.e., In other words, we assume that non-degenerate perturbation theory applies to most realizations of the noise, which requires αV jk E 0 λ,j − E 0 λ,k for most Hamiltonians Ĥλ in the random ensemble Eq. ( 12).Hence, the eigenenergy distribution shows level-repulsion for all state coupled by V .This is for instance satisfied for chaotic quantum systems, which by definition exhibit level repulsion [28][29][30].A multitude of other systems can also be treated within the scope of non-degenerate perturbation theory, such as the hopping of an excitation in a one-dimensional chain with an electric potential ladder [31], or in a molecular network [32].
As we will show with examples in Section V below, if the above conditions are satisfied, the perturbative master equation, Eq. ( 13), captures the ensemble-average dynamics not only on transient time scales, but also on asymptotic time scales (up to systematic errors).This can be understood as a consequence of the interplay between the accumulated errors in the phases from the perturbative expansion, and the ensemble-averaging induced dephasing.On the one hand, the finite perturbative expansion leads to an approximation of the exact phases and eigenvectors of the system.In a closed system these errors in the phases accumulate over time, and eventually any phase relation between the perturbative dynamics and the actual dynamics is lost at sufficiently long times.On the other hand, the ensemble-averaging induced dephasing leads to an overall diminution of the phases and keeps the accumulated error at a finite value.Hence, as long as the dephasing process is fast enough to keep the error from perturbation theory finite, the effective nonunitary dynamics arising from the corresponding perturbative master equation only incur a systematic error.
When the conditions for non-degenerate perturbation theory are not satisfied, the perturbative master equation captures the dynamics at least on short time-scales in the sense discussed in Refs.[4,11,33].Indeed, the first-order in time approximation of Eq. ( 13) yields This equation corresponds to the short-time master equation [4] with Lλ = ( Ĥλ − Ĥ )/( ω 0 ), which captures the effective ensemble-averaged dynamics of any random Hamiltonian ensemble at times t 1/ω 0 , which corresponds to the Gaussian decay regime [? ].

C. Symmetric noise distributions and decoherence rates
Physical noise distributions q ij (∆ jk ) are typically symmetric [34], such as, e.g., Gaussian, Lorentzian or uniform distributions.The non-degenerate constraint E 0 λ,j − E 0 λ,k = (ε j − ε k ) + ω 0 ∆ jk α for most realization of the noise then implies that the distribution must be sufficiently peaked ; if the distribution is too wide, the probability for having degenerate levels could become too high, restricting the validity of the perturbative master equation to short times (Gaussian decay regime), c.f. Section IV B. The decoherence rates Eqs.( 15)-( 17) can then be expressed in terms of the real part (even cumulants) of the characteristic function φ * jk (t) and the central values The zeroth-order function Υ jk separates into an imaginary coherent (energy shift) contribution and a realvalued, time-dependent decoherence rate.The first-order time-dependent decoherence rate Γ jkrj (t) is real-valued, whereas Γ jk (t) is modulated by a fast decaying, complex oscillating function 1−e − it (εj −ε k ) φ * jk (ω 0 t).Note that by the Riemann-Lebesgue lemma [35,36] the characteristic function vanishes at large times, i.e., φ * jk (t) t→∞ −→ 0. We remark that if the distribution q jk is not symmetric, such as, e.g., a general Lévy distribution, one obtains time-dependent energy shifts, i.e., Im [Υ jk ] = Im [Υ jk ] (t) [4].In addition, the energy terms [E 0 j / − E 0 k / ] in the decoherence rates, Eqs. ( 25) and (26), must be replaced by −Im [Υ jk ] (t).

D. Asymptotic state
Given the decoherence rates ( 24)-( 26), the asymptotic state of the perturbative master equation can be expressed in terms of the ensemble-averaged Hamiltonian Ĥ .Indeed, the latter is at large times in the kernel of the perturbative master equation: upon insertion of Ĥ = Ĥ 0 + αV = j (ε j + ω 0 λ j ) |j j| + α V into Eq.( 13), and using the rates ( 24)-( 26), we see that the first commutator vanishes, the second and third sum cancel each other for t 0 because φ * jk (t) → 0, and the fourth sum is of order α 2 (note that the fourth sum vanishes exactly for identically and independently distributed random variables (i.i.d.) because in this case Γ jkrj (t) ≡ 0).The ensemble-average Hamiltonian thus characterizes the asymptotic states up to first order in α.Note that this does not mean that there is a unique asymptotic state; rather, as we will show in the examples below, the perturbative master equation describes asymptotically a dephasing process in the eigenbasis of Ĥ = n E n |n n|.Hence, an initial state ρ0 will asymptotically be projected onto the basis {|n }, As will become clear with the examples in Section V, this means that the coherences (off-diagonal elements ρ jk (0)) of the initial state expressed in the eigenbasis of Ĥ0 , which are coupled by the perturbation V jk , are partially protected from the averaging-induced dephasing.

V. APPLICATIONS
As applications of the perturbative master equation for diagonal noise we consider three examples.(A) First, we study the dynamics of a single qubit which allows for a geometric intuition in terms of Bloch vectors of the dephasing dynamics arising from the ensemble-average.(B) Second, we consider a one-dimensional lattice model with uncorrelated on-site disorder and study the effect of short and long-range interactions on the coherences in the asymptotic state.(C) Third, we investigate the effect of strong correlations on the effective dynamics of the coherences in a many-particle boson model.

A. Geometrical interpretation: Qubit with Gaussian energy disitribution
We begin with two-level systems (d = 2) with static noise in the energy difference, that physically, for instance, may describe an ensemble of spins 1/2 precessing in a static, spatially inhomogeneous magnetic field.It is known that this noise leads to the inhomogeneous broadening of the line-width in spectroscopy experiments, and characterizes the decay of the total magnetization on the timescale T * 2 [9].While the effect of the ensemble-average decoherence can in this case be cancelled out by inverting the spin-precession direction with a spin-echo sequence [9], the latter may also be subject to noise making the inversion incomplete [37].Thus, understanding the exact dynamical role of the static noise is of relevance.
This scenario could be measured in a free-inductiondecay (FID) experiment with nuclear spins [9].After the π/2-pulse used to flip the magnetization in the transverse plane has been applied, the weak radio-frequency transverse field is not turned off, but only turned offresonance.Then, in the rotating frame and considering Gaussian distributed spatial inhomogeneity of the polarization field, we obtain a Hamiltonian ensemble as described in Eq. (28).The FID decay of the magnetization M x on the T * 2 time-scale in the rotating frame then is proportional to Tr [ρσ x ] = Re ρ12 , and should not decay completely due to the presence of the transverse off-resonance field, c.f. Fig. 2. A similar protection of coherences from dipole-dipole interaction T 2 decay by anomalous resonance conditions for the transverse field was discussed in [27,38].
The effective dynamics, Eqs. ( 32), (33), are best understood by expressing the master equation (13) in the eigenbasis of Ĥ by applying the rotation R = e −iθ/2σy to all operators, i.e., Ô → Ôr = R. Ô. R † with θ = arctan[λ 0 /(2α)] the angle between Ĥ and the z-axis.Neglecting the short-time contributions to the rates (γ(t) ≈ −ω 0 σ 2 t/λ 0 ) and the terms proportional to ασ 2 /λ 0 (since |λ 0 | σ + α), we to obtain a pure dephasing equation This differential equation can be solved analytically and yields We obtain immediately from Eqs. ( 36),(37) that the asymptotic state is Applying the inverse transform R † ρr R, we find that the asymptotic state is the projection of the initial state onto the eigenbasis of the ensemble-averaged Hamiltonian Ĥ , c.f., Eq. ( 27).In other words, in the eigenbasis of Ĥ0 the ensemble-averaged dynamics result in a dephasing process towards the basis defined by the eigenvectors of Ĥ .Thus, initial coherences defined with respect to the eigenbasis of Ĥ0 do not decay asymptotically as illustrated in Fig. 3.

B. Effect of long-range couplings: Lattice in the strong bias limit with fully uncorrelated disorder
We consider a one-dimensional tight-binding model with d sites with on-site disorder and a constant potential ladder as illustrated in Fig. 4. The random Hamiltonian ensemble is parametrized as where λ j are identically independently distributed (i.i.d.) Gaussian (c.f.Eq. ( 29)) dimensionless random variables with x = 0, 1, 3.This model for instance describes s-band electron transport in the presence of a strong electrical potential [31], and corresponds to an Anderson-type model [16] in the strong bias limit.Experimentally, it can be simulated using atom-optics simulators [39,40] or photonic wave-guides [41][42][43].
The eigenvalues of the unperturbed part Ĥ0 λ = Ĥs λ + T are E 0 j,λ = ω 0 (jT + λ j ) (i.e.ε j = ω 0 jT ).As was derived in Section III, the corresponding first-order in α perturbative master equation, Eq. ( 13), is fully characterized by the characteristic function, Eq. ( 20), Using the latter and Eqs. ( 24)-( 26) we obtain for the decoherence rates where we neglected the fast decaying contribution to γ jk (t).Here Γ jkrj (t) = 0 because the eigenvalues are i.i.d.Note the factor 2 in Eqs. ( 42) and ( 43) as compared to the single qubit case, Eqs. ( 30) and (31), which is due to the i.i.d.condition.As a numerical example, we consider a system of d = 30 sites with σ = α and T = 10α, and a broad, centered Gaussian initial state with N the normalization constant and p Ga a Gaussian distribution with average (d+1)/2 and variance √ d.This could for instance model a photo-excitation in a correlated thin-film transition metal-oxide heterostructure with a strong intrinsic electrical field [31].
In Fig. 5, we show the evolution of the total coherence for VNN , Vx (x = 0, 1, 3), and V = 0.The dynamics are obtained by numerical integration of the perturbative master equation ( 13) with the 'NDsolve' routine in Mathematica 11.0.At short times, the coherences decay, independently of the coupling potential, as ∼ e −σ 2 ω 2 0 t 2 under the action of the rate Υ jk .At large times, for any non-vanishing potential V , the decay eventually stops, and a finite amount of coherence remains in the asymptotic state, which is in agreement with Eq. ( 27).For long-range potentials the total coherence of the asymptotic state, c ∞ = lim t→∞ c(t), is larger than for shortrange potentials.Indeed, for a fully connected network (x = 0) the total coherence in the asymptotic state is c ∞ = 15.1 × 10 −3 , for slowly decreasing coupling (x = 1) we obtain c ∞ = 4.1 × 10 −3 , and for dipoletype coupling (x = 3) we have c ∞ = 1.6 × 10 −3 .Interestingly, the nearest-neighbours coupling converges to c ∞ = 1.4×10 −3 , which is close to the value for the dipole coupling.
As proof of consistency of the perturbative master equation approach, we computed the full dynamics by numerical exact averaging using the numerical solver from the python Qutip package 4.3.0[44].For the long-range interactions x = 1 and x = 0, we found that 10 6 realizations are sufficient for the convergence of the numerical averaging.The results of the master equation and direct averaging well agree as shown in Fig. 5.For the shorter range interactions (x = 3, VNN and V = 0), 10 6 realizations were insufficient.In Fig. 5 we show the deviations for V = 0 for which we know that the master equation is exact [4].Note that on the same computer the numerical integration of the master equation was ten times faster than the direct numerical averaging with 10 6 realizations.As further test we verified that the fidelity between the density matrix from the numerical computations and the master equation is larger than 0.999 at all times.Furthermore, we found that the relative purity Tr ρ 2 /Tr ρ 2 Num shows deviation of up to 6% in the time range where the decay enters the asymptotic state (1.5 t 2.5), but eventually converges to the same value (for x = 0, 1).This deviation occurs due to the approximation in the time-dependence of the decoherence rates Eqs.( 24)- (26).46), for different dipole-type coupling potentials with exponent x = 0, 1, 3 (c.f.Eq. ( 41)), nearest-neighbour coupling (NN) and no coupling (V = 0).The longer-range the coupling, the more coherences are present in the asymptotic state characterized by Eq. ( 27) (horizontal gray dashed lines).Numerical integration of the perturbative master equation, Eq. ( 13) with rates Eq. ( 42)-( 44), (solid lines) agrees with the direct numerical averaging (dots) with 10 6 realizations of the disorder for x = 0, x = 1.For x = 3, NN and V = 0 (only V = 0 is shown for visual purposes) the asymptotic state is not reached due the statistical error ∼ √ 10 6 .The parameters are d = 30, σ = α, T = 10α and the initial Gaussian state, Eq. ( 45), has width √ 30 and average value 31/2.
Overall we find that the dephasing induced by diagonal disorder does not, in the presence of couplings V between the eigenstates, lead to a full decay of the coherences.This is in stark contrast to dynamical diagonal noise such as considered in the Hacken-Strobl model [45] or for homogeneous broadening descriptions [46].Moreover, the coherences are most efficiently protected from the disorder-induced decoherence by long-range couplings.This can be understood by analogy to the qubit case studied in Section V A where we demonstrated that the ensemble-averaging leads to an effective dephasing in the eigenbasis {|n } of the ensemble-averaged Hamiltonian Ĥ (c.f.Fig. 3).For V = 0, the basis {|n } is equal to the quantization basis {|j }, and we thus have a pure dephasing process so that all coherences vanish.However, the stronger the potential V , the more the eigenstates of Ĥ will deviate from |j , and thus the more the dephasing basis {|n } differs from the quantization basis.Consequently, a larger amount of the coherences of the initial state do not decay asymptotically.
This result could be of relevance for a better understanding of the interplay between disorder and quantum coherent transport (which relies on the coherences) in strongly connected networks such as the Fenna-Matthew-Olson photo-synthesis molecular complex [47], assemblies of ultra-cold Rydberg atoms [48], strongly-correlated materials [31,49], superconducting circuits [50] or photonic circuits [51]. .

C. Many-body (bosons) dynamics with strongly correlated noise
We consider an asymmetric double-well potential with N interacting bosons (see illustration Fig. 6) described in terms of a tilted Bose-Hubbard model [52][53][54]).Such a setting can be experimentally implemented with ultracold bosons in optical lattices [55,56] and was studied, e.g., in the context of quantum Chaos [57], superfluidity [58], and distinguishability [59].Static diagonal noise can either arise from random on-site interaction and/or from fluctuations in the tilt between the left and right wells.In an experiment with optical lattices, the former may arise from imprecisions in the Feshbach resonances used to fix the magnitude of the interaction [60,61], while the former comes from fluctuations of the trapping potential.Considering one or the other is equivalent from a formal point of view, and we arbitrarily choose the interaction to be random.
Firstly, we fix the notation and write the Bose-Hubbard Hamiltonian in the absence of noise as with the unperturbed kinetic term + ω 0 U 0 2 NL ( NL − 1) + NR ( NR − 1) , and the perturbation Here âL , âR , â † L , â † R are the bosonic annihilation and creation operators of the left (L) and right (R) wells respectively, and the number operators NL = â † L âL , NR = â † R âR .The tunnelling rate is denoted as ω 0 J ∈ R + , the on-site interaction reads ω 0 U ∈ R and the tilt between the two wells is given by ω 0 T ∈ R. The eigenvalues of the unperturbed Hamiltonian Ĥ0 can be writ-ten as E 0 m = ω 0 (β m U 0 + χ m T ) with the constant factors β m = 1/2m(m − 1) + 1/2(N − m)(N − m − 1) and χ m = (2m − N ), m = 0, . . ., N fixed by the number of bosons N in the system and using the Fock basis ordering |N, 0 , |N − 1, 1 , . . .|0, N .For example for N = 3 we have Adding a random noise δU to the interaction, we replace U 0 → U 0 + δU .Consequently, the eigenvalues of the unperturbed random Hamiltonians are given by E 0 m,λ = ω 0 (β m U 0 + χ m T ) + ω 0 λ m , with λ m = δU β m .The noisy eigenenergies are thus strongly correlated as they are characterized by a single random variable δU , i.e., the eigenvalues are not independently distributed, as opposed to the previously studied model of a onedimensional potential ladder with on-site disorder.
Assuming a Gaussian distribution p Ga (δU ) of mean value 0 and variance σ, the characteristic function, Eq. ( 20), evaluates to φ * jk (ω , and the decoherence rates Eqs.( 24) -( 26) read Hence, second order processes described by Γ jkrj (t) now play a role because the eigenvalues are strongly correlated.To guarantee that the eigenvalues E 0 m,λ of the unperturbed random Hamiltonian remain non-degenerate for most realizations of the noise, we require |T | σ.For the perturbation to remain small, we further impose the condition αJ < |U 0 + T |.
Hence, in addition to coherences being present in the asymptotic state due to the coupling V , the symmetry of the Hamiltonian gives rise to slowly decaying coherences.This effect could be exploited to generate longlived coherences of many-body states in systems subject to generic on-site noise.

VI. CONCLUSIONS
We derived general perturbative master equations to describe the ensemble-averaged dynamics of quantum systems that are subject to static noise (which includes disorder, the coupling with a classical environment, or slowly drifting experimental parameters).The pertubative master equations not only provide a description of the ensemble-averaged dynamics in terms of physically interpretable operators, energy shits, and (timedependent) decoherence rates, but also require much less computational resources to numerically integrate as compared to direct numerical averaging of the dynamics.
We treated in details systems described by Hamiltonians with diagonal static noise, and a noise-free perturbative coupling potential.In the range of parameters where perturbation theory converges for most realizations of the noise, the first-order master equation describes the ensemble-averaged dynamics on all time scales up to systematic errors.The effect of the ensemble-averaging can be understood as a dephasing (decoherence) process in the eigenbasis of the average Hamiltonian, which was illustrated in the Bloch sphere for a two-level system with a random energy splitting.Thus, the asymptotic state is the projection of the initial state on the eigenbasis of the latter.This was numerically verified using the perturbative master equation and direct numerical averaging for several examples.
The generality of our approach was shown with a onedimensional tight-binding model with on-site disorder and a Bose-Hubbard double-well Hamiltonian with onsite interaction noise.In the first example, we showed that the longer range the perturbative coupling potential is, the more coherences are protected from the dephasing and remain in the asymptotic state.In the second example, we showed that strong correlations in the noise distribution result in a slow decay of the coherences, which are thus partially protected from the averaging-induced decoherence.For all scenarios we discussed various experimental setups to measure these effects.Our work suggest that the later could be exploited to protect coherences for quantum applications.Conversely, static random distribution could be introduced on purpose to give rise to a desired type of dephasing, and, for instance, drive a quantum system into a target non-equilibrium quantum state.
and the master equation matrix, Eq. ( 9), then reads Appendix B: First-order Perturbation equation Here we show the computations for deriving the master equation (13).For more details see [11].
a. Zeroth order α 0 : The zeroth order dynamical matrix captures the dynamics of solely the unperturbed part, and describes the evolution of the coherences in the eigenbasis {|j }.We recall that Û 0 λ (t) = e − it Ĥ0 λ = j e − it E λ,j |j j| and derive the ensemble-average value where the ensemble-averaged phase factors and thus Note that one must remain careful here because even tough ϕ −1 jk (t) = ϕ * jk (t), in general ϕ −1 jk (t) = ϕ * jk (t).b.First order α 1 : From Eqs. ( 6) and (8) where we defined for simplicity of reading and where for f jkrj the ordering of the indices refers to the ordering of the indices of the phase factors.Note that in the above above, we separated out the terms for which the phase factors are equal to the identity prior to the ensemble-averaging because in general ϕ jk (t) δ j,k = ϕ jj (t) = 1.
Using the zeroth Eqs.(B1) and (B3) and the first order Eq.(B4) for the dynamical matrix, we can compute the first-order expansion of Q , Eq. ( 9), which fully characterizes the perturbative master equation.To zeroth order we obtain The first order yields The symmetry Γjkrj = − Γ * kjjr is obtained by permuting the first two indices, jk → kj, and the last two indices, rj → jr, and reflects the Hermiticity of the density matrix.
In order to derive a master equation in Lindblad form following the general method described in [4], we should now expand Q 1 (t) in a basis of Hermitian traceless operators and collect the terms for the coherent and incoherent parts.However, this proves to be, in general, technically rather involved.We prefer to go back to the Hilbert space representation using the identity ρ(t) = Let us derive the decoherence rates γ jk (t),Eq.( 16), and Γ jkrj (t), Eq. ( 17), in the case that the joint-noisedistribution q jk (λ j − λ k ) for any pair of variables λ j , λ k has vanishing odd cumulants (except the central value), or in other words, the probability density functions q jk are symmetric around their central value.In this case, considering furthermore that the probability distributions shall be continuous and smooth, we know that the associated characteristic function vanishes in the limit of large times, φ jk (t) = dλ j λ k q jk (λ j − λ k )e −i(λj −λ k )t t→∞ −→ 0.

Figure 1 .
Figure 1.Illustration: Static noise in the on-site energies of a chain leads to a decay of coherences ρ jk (off-diagonal elements) of the ensemble-averaged density matrix, but longrange coupling V can prevent a complete decay in the asymptotic state ρ∞ .

Figure 3 .
Figure 3. Illustration in the Bloch sphere of the ensembleaveragind induced dephasing process (blue line) Eqs.(32)(33) from static diagonal noise.The asymptotic state is the projection of the initial state onto the average Hamiltonian Ĥ (green line) and has non-vanishing coherences (x and y components).For visualization purposes the initial state is different from Fig. 2.

Figure 4 . 1 j=1
Figure 4. Illustration of a one-dimensional potential ladder with short and long-range coupling V , a potential step T and on-site disorder λj (c.f.random Hamiltonians defined in Eq. (39)).

Figure 6 .
Figure 6.Double-well potential with a tilt T between the wells left (L) and right (R) filled with bosons (blue dots) with onside interaction U and tunnelling barrier J as an illustration of the Bose-Hubbard Hamiltonian(47)..