Breakdown of quantum-to-classical correspondence for diffusion in high temperature thermal environment

We re-consider the old problem of Brownian motion in homogeneous high-temperature thermal environment. The classical theory implies that the diffusion coefficient does not depend on whether the thermal fluctuations are correlated in space or disordered. We show that the corresponding quantum analysis exhibits a remarkable breakdown of quantum-to-classical correspondence. Explicit results are found for a tight binding model, within the framework of an Ohmic master equation, where we distinguish between on-site and on-bond dissipators.

High temperature classical Brownian motion is described by the Langevin equatioṅ where f is white noise of intensity ν and, and η = ν/(2T ) is the friction coefficient. Consequently the diffusion coefficient in space is D = T /η. The white noise arises in general from a fluctuating potential, namely f = −∂ x U(x, t), but the classical result does not depend on its spatial correlation scale . In fact the common practice is to assume =∞, aka the Caldeira-Leggett model [1,2], meaning that f is independent of x. But in the quantum treatment does show up in the analysis, because it determines the lineshape of the stochastic kernel W(k|k ) for scattering from momentum k to momentum k. Namely, the width of the kernel (∼2π/ ) has implication on the transient decoherence process [3][4][5]. Yet, one does not expect that this lineshape will have any effect on the long time spreading. The argument is simple: on the basis of the central limit theorem successive convolutions should lead to a result that does not depend on the -dependent lineshape of the stochastic kernel, but only on its second moment, which is characterized by ν. Consequently robust quantum-to-classical correspondence (QCC) is expected at high temperatures. Such QCC can be regarded as an implication of the Thomas-Reiche-Kuhn-Bethe-Wang sum rule [6], or as an extension of the restricted QCC principle [7,8].
Scope and outline.-In the present work we show that independence of D is a fallacy, and consequently D acquires a non-universal dependence on the temperature. This is demonstrated for a Brownian motion on a tight binding chain as in [9][10][11][12][13][14][15][16][17][18][19][20]. The interest in such models has been further motivated by the puzzle of excitation transport in photosynthetic light-harvesting complexes [21][22][23][24][25][26][27][28][29][30], and it is somewhat related to past studies of motion in washboard potential [31][32][33]. We derive the dynamics from an Ohmic master equation, and further illuminate the results using an effective rate equation. We emphasize that our results have nothing to do with the well known studies of quantum dissipation in low-temperatures [34][35][36][37], where non-classical effects are related to the failure of the Markovian approximation (Ohmic fluctuations in low temperatures become strongly correlated in time). See [SM] for a short overview. Before going into details we present the model system and the main results.
Model and main results.-The isolated chain is defined by the Hamiltonian H (c) = −c cos(p) that describes a particle or an exiton [38][39][40] that can hop along a onedimensional chain whose sites are labelled by an integer index x. The operators e ∓ip generate one-site displacements. The inverse hopping frequency 1/c corresponds to the effective mass of the particle. The dynamics is governed by an Ohmic master equation for the probability matrix Following [41] we distinguish between 3 types of dissipators: (a) The Caldeira-Leggett X-dissipator L (X) that corresponds to non-disordered ( =∞) bath. (b) The Sdissipator L (S) due to uncorrelated noisy sites. (c) The B-dissipator L (B) due to uncorrelated noisy bonds. With the X-dissipator we find [we shall add a reference if it turns out that this result is already known]: where I n is the modified Bessel function. This result is exact to the extent that the (Markovian) Ohmic master equation can be trusted. For low temperatures (in the sense T c) one recovers the standard result D = T /η which applies for non-relativistic p 2 dispersion. For high temperatures, the diffusion is composed of a term that originates from the coherent hopping between the sites ( ) and a term that arises due to non-coherent hopping, induced by the bath (⊥). Accordingly a general expression for the diffusion is given by: for X/S-coupling C ⊥ =0. Similarly, by convention, we have for X-coupling C =1. With the same conventions we get for B-coupling C =1/6, and for S-coupling C =1/2. In contrast the A coefficients are not a matter of convention. Rather they reflect the thermalization and the spreading mechanism, and hence indicate quantum manifestation: The dependence of D on the temperature is plotted in Fig. 1 and will be further analyzed below. For sake of comparison we plot also the naive expectation D ∝ v 2 , with v = c sin(p), where the average is over the canonical distribution. The high-temperature dependence is We shall explain that with S/B dissipators the classical evaluation of D leads to wrong results, while the quantum version leads to Eq. (4). Nevertheless, unlike the X-coupling case, for local (S/B) dissipators the steady state ρ is canonical only in second order, meaning that Lρ ∼ O(β 3 ). This is a general result of the steady state in the Ohmic approximation [42]. If we ad-hock correct the transition rates to get agreement with Boltzmann, the results for the A-s are modified as follows [SM]: We emphasize again that the value of A is a sensitive probe that is affected by the line-shape of the spreading kernel. Therefore its precise value is non-universal but depends on the weights of the quantum transitions.
The Ohmic dissipator.-The isolated chain is defined by the H (c) Hamiltonian. The X-dissipation scheme involves a single bath, with interaction term −W F , where W is the position operator x, and F is a bath operator that induces Ohmic fluctuations with intensity ν. More generally we assume a disordered thermal environment that is composed of numerous uncorrelated baths such that the interaction term is α W α F α , where α labels different locations. For Sdissipation W α = |x α x α |, leading to a fluctuating potential that dephases the different sites. For B-dissipation W α = |x α +1 x α | + h.c., which induces incoherent hopping between neighbouring sites. The Ohmic dissipator L (X/S/B) ρ takes the form where η = ν/(2T ) is the friction coefficient, and The friction terms represent the response of the bath to the rate of change of the W α . For X-dissipation V = c sin(p) is the velocity operator. If we treat the the friction term of Eq. (9) in a semi-classical way, the expression for the dissipator in the Wigner phase-space representation ρ w (R, P ) takes the familiar Fokker-Plank (FP) form with v = c sin(P ), namely, which is a sum of momentum-diffusion and momentumdrift terms. For the sake of later reference we have added to the friction force (−ηv) a constant field f 0 . The X-dissipator leads to canonical steady state (SS) for any friction and for any temperature. This is not the case for S/B-dissipation, for which the agreement of the SS with the canonical result is guaranteed only to second order in η. The reason for that is related to the proper identification of the "small parameter" that controls the deviation from canonical thermalization. The X-dissipator induces transitions between neighboring momenta, and therefore the small parameter is ∆/T , where the level spacing ∆ goes to zero in the L → ∞ limit, where L is the length of the chain. But for local baths, the coupling is to local scatterers, that create transitions to all the levels within the band. Therefore the small parameter is c/T , and canonical thermalization is expected only for c/T 1.
Semiclassical analysis for for X-dissipation.-We shall argue below that for X-dissipation the classical dynamics that is generated by L FP is exact for the purposed of the A coefficient evaluation. The calculation of the velocity-velocity correlation function v(t)v(0) requires a rather complicated recursive procedure [SM]. The final step is to get the diffusion coefficient via Technically there is a shortcut that can be used in order to get D without going through the complicated calculation of the correlation function. The way to go is to add in Eq.(11) a constant field f 0 . To find the SS is rather simple. For weak field one obtains v = µf 0 , where µ is the mobility. Then the diffusion coefficient is deduced from the Einstein relation, namely, D = µT . The same result Eq. (3) is obtained in both methods [SM]. Semiclassical analysis for S-dissipation.-In the classical treatment x is regarded as as a continuous coordinate, and therefore we write W α = u α (x) = u(x−x α ), that involves a short-range interaction potential u(r). Let us denote by ν (S) the fluctuations of the potential These fluctuations have the same intensity at any x because we assume that the x α are homogeneously distributed. Then it follows that the fluctuations of the force, that is derived from the stochastic potential, have an intensity ν = (1/ ) 2 ν (S) , where is the correlation scale. We set ∼ a where a=1 is the lattice constant. Then it follows automatically that η = ν/2T is still the friction coefficient, as in the case of X-dissipation. See [4] for details. So in the classical description we get the same Langevine equation, irrespective of the correlation distance that is determined by the width of u(r).
Semiclassical analysis for B-dissipation.-Using the same prescription, and ignoring commutation issues, we get W α = [2 cos (p)]u α (x). This means that motion with momentum |p| ∼ π/2 is not affected by the baths. We shall see later that this is an artifact of the classical treatment, and does not hold for the quantum dynamics. Still, the classical perspective provides some insight that helps to clarify how the second term in Eq.(4) comes out. The equations of motion that are derived from the full Hamiltonian are of Langevin-type: For infinite temperature the F α are uncorrelated white noise terms, with some intensity proportional to ν (B) . Therefore we get from Eq.(14) diffusion in p with coefficient where ≈ a and a=1. The latter term, after momentum averaging, is responsible for getting the D ⊥ term in Eq. (4). For a particle that moves with constant momentum p, ignoring the variation in p, the velocityvelocity correlation decays as exp(−ν x t) due to this xdiffusion. This leads to an extra Drude term D = v 2 /ν x that diverges at p=π/2. However, taking the variation of the momentum into account, this divergence has zero measure, and the final result [SM] is finite, leading to the first term in Eq. (4) with C = 0.49. For finite temperature the fluctuations gain non-zero Quantum dynamics.-The quantum evolution is generated by L of Eq.(2) with the dissipators of Eq.(9), and it can be written as sum of Hamiltonian, noise and friction terms, namely L = cL (c) + νL (ν) + ηcL (η) . The elements of the super-vector ρ are given in the standard representation by ρ(R, r) ≡ R + r/2|ρ|R − r/2 , and in Dirac notation we write ρ = R,r ρ(R, r) |R, r . The super-matrix L is invariant under R-translations, and therefore it is convenient to switch to a Bloch representation ρ(q; r) where L decomposes into q blocks. In the q subspace we have the following expressions [SM]: The subscripts X/S/B distinguish the different coupling schemes, and D ⊥ = |r+1 r| is the displacement operator in r space.
The quantum analysis.-To obtain the diffusion coefficient, we consider the spectrum of L for a finite system of L sites. In the Bloch representation the equation Lρ = −λρ decomposes into q-blocks. For a given q we have a tight binding equation in the |r basis. For example L (c) induces near-neighbor hopping in r. The eigenvalues for a given q are labeled λ q,s , where s is a band index. The long-time dynamics is determined by the slow (s=0) modes. Specifically, the diffusion coefficient is determined by the small q expansion λ q,0 ≈ Dq 2 .
The NESS eigenvector belongs to the q=0 block, and for η=0 it is given by |r=0 . Non-zero q and η can be treated as a perturbation. The key observation is that in order to get an exact result for D it is enough to use second-order perturbation theory in q. The outcome of this procedure is the analytical expression Eq.(4) with the associated results for the A coefficients. See [SM] for extra technical details.
Wigner phase-space picture.-The propagation of the Wigner distribution function ρ w (R, P ) is generated by a kernel L(R, P |R 0 , P 0 ) that is obtained from Eq. (15) in a straightforward manner via Fourier transform [SM]. For simulations of the long time spreading it is enough to approximate L in a way that is consistent with second order perturbation theory in q. As explained in the previous paragraph, such approximation provides an exact result as far as D calculation is concerned. Replacing sin(q/2) by (q/2), the L (c) term by itself generates classical motion in the X direction with velocity v = c sin(P ). In the quantum calculation this motion is decorated by a Bessel function, but D is not affected. The cos(q) in the L (ν B ) , after expansion to second order and Fourier transform, leads to an x-diffusion term, that is responsible to for the C ⊥ contribution in Eq.(4). As far as this term is concerned, there is no difference between the quantum and the classical picture, and therefore we ignore it in the subsequent analysis. The cosine factors in the other dissipators can be replaced by unity. The reason is as follows: by themselves those cosine terms do not lead to any diffusion; only when combined with the L (c) term they lead to the Drude-type C contribution in Eq. (4); the L (c) is already first order in q; hence no need to expand the cosines beyond zero order.
Effective rate equation.-With the approximations that were discussed in the previous paragraph (excluding for presentation purpose the trivial R diffusion), we find that the evolution of the Wigner function is generated by a stochastic-like kernel L(R, P |R 0 , P 0 ) = W(P |P 0 )δ(R − R 0 ). The explicit expressions for infinite temperature (η=0) are: These are the transition rates (P = P 0 ), while the diagonal elements of W are implied by conservation of probability. For X-dissipation Eq.(16) describes local spreading of momentum which is in complete correspondence with the classical analysis. For S-dissipation Eq.(17) describes quantum diffractive spreading. In the latter case, if the dynamics were treated classically one would obtain the same result as for X-dissipation, namely Eq. (16), with prefactor of order unity that can be by re-scaled to unity by adopting the appropriate convention for the definition of ν. In other words: the coupling strength to the bath should be re-defined such that ν is the second-moment of W(P |P 0 ) irrespective of the lineshape. Similarly, if the dynamics were treated classically for the B-coupling, one would obtain Eq.(16) multiplied by 4 cos 2 (P ), as implied by the classical analysis. The result for W for finite temperature, in leading order in η (which serves here as a dimensionless version of the inverse temperature) can be written as where E(P ) = −c cos(P ). More precisely, if we incorporate the L (η) term of the Ohmic master equation, we get Eq. (19) with e x → (1 + x). This reflects the well known observation that the Ohmic approximation satisfies detailed balance to second order in η. Accordingly the Ohmic SS agree to second order with the canonical SS, namely, ρ SS (P ) ∝ exp[−E(P )/T ]. Stochastic simulations.-In Fig.1 we test the analytical approximation Eq.(4) against exact numerical calculation that is based on the effective rate equation. The diffusion coefficient D is calculated using Eq. (12). The momentum spreading kernel is K(t) ≡ exp(Wt), and the velocity is v P = c sin(P ). Accordingly In the tight-binding framework we have the identification m → 1/(ca 2 ), where a is the lattice constant. There is a crossover to standard QBM as θ ≡ T /c is lowered. It is illuminating to summarize this crossover in terms of mobility. Using the Einstein relation we get The standard result is the first term with B(θ) = 1, while Eq.(4) implies that B(θ) ∝ [(1/θ) 2 +A (1/θ) 4 ] for large θ.
Additionally the C(θ) term is due to bath-induced incoherent hopping. The A coefficients that control the temperature dependence of B(θ) and C(θ) provide a way to probe the underlying mechanism of dissipation, and to identify the high-temperature fingerprints of quantum mechanics.

====== [1] Overview
There is a vast literature regarding Quantum Brownian Motion (QBM). The prototype so-called Caldeira-Leggett model corresponds to the standard Langevin equation where the dispersion relation is v = (1/m)p. In the tight-binding framework we have the identification m → 1/(ca 2 ), where a is the lattice constant. The standard QBM model features a single dimensionless parameter, the scaled inverse temperature β, which is the ratio between the thermal time 1/T and damping time m/η. In the lattice problem one can define two dimensionless parameters Accordingly β = α/θ. In our model we set the units such that a = 1, hence, disregarding 2π factor, our scaled friction parameter η is the same as α. The regime diagram of the problem is displayed in Fig.S1, and further discussed below. It contains both CBM-like regime, where memory effects are either not expressed or appear as a transient, and QBM regimes where the dynamics is drastically different. The standard analysis of QBM [35] reveals that quantum-implied memory effects are expressed in the regime β 1, where a transient log(t) spreading is observed in the absence of bias, followed by diffusion. Most of the following quantum dissipation literature, regarding the two-site spin-boson model [43] and regarding multi-site chains [12,33], is focused in this low temperature regime where a transition from CBM-like behavior to over-damped or localized behavior is observed, notably for large α of order unity.
Our interest is focused in the α, β 1 regime. This regime is roughly divided into two regions by the line θ ∼ 1, see Fig.S1. Along this line the thermal de-Broglie wavelength of the particle is of order of the lattice constant, hence it bears formal analogy to the analysis of QBM in cosine potential [32], where it marks the border to the regime where activation mechanism comes into action. In our tight binding model we have a single band, hence transport via thermal activation is not possible. Rather, in the θ > 1 regime the momentum distribution within the band is roughly flat, We distinguish between the Classical-like Brownian Motion (CBM) region; the low-temperature QBM region where memory effects dominates; and the high-temperature QBM region that has been discussed in this work. Note that below the dashed diagonal line β = 1 memory effects should be taken into account. There is a range of temperatures in the QBM region where they manifest as a transient [35]. (b) The scaled mobility µ/µ0 where µ0 = 1/η versus θ, based on the analytical results that have been obtained for diffusion in the X/S coupling schemes. The result is independent of η. We also add the result for the B coupling scheme that approaches the finite asymptotic value µ∞ = 2η (horizontal line). In the latter case η = 0.3 has been assumed. Note that the S/B analytical results are applicible only in the θ > 1 regime.

====== [2] Classical Brownian motion on a lattice
We consider classical equation of motion for a Brownian particle that has dispersion as in a tight-binding chain, with a coupling to a non-disordered fluctuating field. A fully-quantum version of this system was studied by [12,33], with focus on low temperature QBM regime. In this section we find the steady state in the regime where the Ohmic master equation is valid. Setting the lattice constant to be unity, the Hamiltonian is: The steady state for p is solved by inserting Eq.(S-3) to Eq.(S-4). Changing notation p → ϕ, and u(ϕ) = f 0 − ηc sin(ϕ), and D ϕ = (1/2)ν, one get the equationφ = u(ϕ) + f (t), with the associated Fokker-Planck equation This equation describes motion in a tilted potential The steady state solution is where the integration constant C is determined by the periodic boundary conditions ρ 0 (0) = ρ 0 (2π), namely, Simplifying, the final expression for the NESS is Averaging over Eq.(S-4), and using ṗ = 2πI, one obtains where µ is the so-called linear mobility. This result for µ is consistent with direct calculation of D in accordance with the Einstein relation, namely µ = D/T . The direct calculation of D is performed for zero field. For a particle that starts at x=0, the variance of x after time t is: Therefore D = c 2 S 1 , where S 1 is the area of the sine correlation function, as defined in the next section, where we outline its calculation, leading to Eq.(3) in the main text.

====== [3] The sine correlation function
First we recall that for zero field the steady state is an equilibrium canonical state ρ(ϕ) ∝ exp[−W (ϕ)], where W (ϕ) = z cos(ϕ), and z = (c/T ). At equilibrium we have w n ≡ cos (nϕ) = I n (z) I 0 (z) (S-14) We define S n as the area of the sine-sine correlation function s n (t) = sin(n ϕ t ) sin(ϕ 0 ) , namely, Eventually we are interested only in S 1 , but for the derivation we define a full set of sine correlation functions. Explicitly these are written as The average without subscript assumes equilibrium state, while the average with subscript "0" indicates initial condition ϕ 0 and assumes a Langevin picture. The subscript "t" indicates expectation value after time t within the framework of the associated Fokker-Planck picture. Initially we have s n (0) = sin (nϕ) sin (ϕ) = 1 2 In order to find s n (t) at later times, we realize that the it satisfies the same equation of motion as that of sin(nϕ t ) 0 , where 0 indicates any initial state. This is known as the "regression theorem". The adjoint equation for any observable Where we used the completeness relation

====== [4] The Wigner phase space representation
Here we treat (x, p) as extended continuous coordinates and derive the standard Wigner representation for the quantum propagation in the absence of dissipators. The elements of the ρ are given in the standard space representation by ρ x ,x ≡ x |ρ|x . We define r = x − x and R = (x + x )/2, and use super-vector Dirac notations, namely ρ = R,r ρ(R, r) |R, r . The space representation is ρ(R, r) ≡ ρ x ,x , the momentum representation ρ(q, P ) is related by double Fourier transforms, and the intermediate representations are those of Wigner ρ w (R, P ) and Bloch ρ(q; r). For the unitary evolution with U = exp[ict cos(p)], the propagator of the Wigner function in momentum representation is -24) Note that this kernel is properly normalized with respect to the integration measure dRdP/(2π).
With sin(q/2) → (q/2) we get the classical result But quantum mechanically we get In the above sum n runs formally over all the integer and half-integer values. Note that Wigner function on a lattice has support on both integer and half integer lattice points (weight on half integer lattice points is the fingerprint of interference due to superposition of integer lattice locations).

====== [5] The Bloch representation
For an infinite chain the conventional way to define the Bloch representation is to perform R → q Fourier transform of ρ(R, r) for a given r to obtain ρ(q; r). Note that R runs over integer values for r = 0, 2, 4, ... and over half integer values for r = 1, 3, 5.... This definition has a problem if we consider a finite chain with periodic boundary conditions. Still it can be justified after a short transient if L is large enough because distant points in space loose phase correlation (if there was to begin with). For a small ring (small L) this might not be the case. Therefore in a previous work [41] we have defined ad-hock the Bloch representation ρ q (r) as the Fourier transform of x|ρ|x + r . The ad-hock definition differs by gauge transformation (and non-intentionally also by sign) from the conventional definition, and allows to handle correctly the periodicity in both coordinates, namely, also in r. For a small chain, or for a complete investigation of the eigenvalues problem, these phases are important. See for example [5].
Our system is invariant under translations, therefore it is natural to perform the diagonalization of L is the Bloch representation. In practice one can obtain the expressions in Eq.(15) by inspection. As an example let us see how the expression for L (c) is obtained. It originates from i[cos(p), ρ]. In the standard representation its matrix elements are Recall that cos(p) is the sum of displacement operators e ∓ip . In super-vector notations the above expression can be written in terms of operators e ∓i(1/2)q and e ∓iP that induce translations in R and in r respectively. Namely, Thus we can write In the Bloch (q, r) representation this super-operator becomes block diagonal in q.

====== [6] From Bloch to Wigner
In the main text we present in Eq. (15) the Bloch representation L (q) of the dissipators. The transformation to the Wigner representation is essentially a Fourier transform: Note that inner integral transforms r L (q) r 0 to the momentum representation P L (q) P 0 . Note also that W(P |P 0 ) = P L (q=0) P 0 are the Fermi Golden Rule (FGR) transition rates between momentum eigenstates. Commonly the FGR is considered as an approximation, while we have rigorously established that W(P |P 0 ) can be used within an effective rate equation in order to evaluate the exact quantum result for D.
In the main text we use a discrete momentum notation, such that r|P = L −1/2 exp(iP R), etc. Consequently, in the discrete version of Eq.(S-30), the integrand of the drdr 0 integral contains an extra 1/L factor. On the other hand for summation P over momenta the measure becomes [L/(2π)]dP .

====== [7] The Boltzmann dissipators
It is convenient to handle the calculations of the spectrum on equal footing for all the coupling schemes, for both Ohmic and Boltzmann versions of the dissipators. Consequently, in the latter case we have to transform Eq. (19) back from the Wigner representation to the Bloch representation. Namely, r L (q) r 0 = r, q|L|r 0 , q = 1 L P,P0 W(P |P 0 )e iP r−iP0r0 (S-31) Note that this expression does not depend on q, reflecting the δ(R − R 0 ) of the transitions. Making the distinction between the diagonal terms (P = P 0 ) and the off-diagonal terms (P = P 0 ), taking into account that by definition the kernel conserves probability ( P W(P |P 0 ) = 0), one can write In the latter expression it is implicit that P = P 0 has measure zero, so it reflects the contribution of the P = P 0 terms in the discrete sum of Eq.(S-31). For the S and B schemes one obtains W (S) (r, r 0 ) = νI r (z/2)I r0 (−z/2) (S-34) W (B) (r, r 0 ) = 2νI r (z/2)I r0 (−z/2) + ν [I r+1 (z/2)I r0−1 (−z/2) + I r−1 (z/2)I r0+1 (−z/2)] (S- 35) where z = c/T . The high-temperature calculation of the spectrum, Eq.(S-33) is Taylor expanded in z = (c/T ) up to second order using The first order result for L (q) is a q = 0 version of the S/B dissipators that were presented in the main text Eq. (15). Both schemes acquire second-order terms −(3/16)(c/T ) 2 ν |±1 ±1| that are required for the calculation of D, see next SM section. For the S coupling scheme one finds additional second-order terms that are needed for the calculations, namely, −(3/32)(c/T ) 2 ν |1 −1| and −(3/32)(c/T ) 2 ν |−1 1|.

====== [8] Perturbation theory
We use perturbation theory to find the eigenvalue λ q,0 of L (q) , from which we can obtain D as explained in the main text. We regard the Bloch quasimomentum q and the friction η as the perturbation. For q = η = 0 the state |r = 0 is an exact eigenstate that is associated with the eigenvalue λ = 0. Due to the perturbation it is mixed with neighboring |r states. We outline below how we get analytical expressions for λ q,0 to any order in q and η. In practice we go up to second order.
In the following we demonstrate how we perform perturbation theory for the X-coupling scheme. The same method is used for the S/B coupling schemes either with the Ohmic dissipators or with the Boltzmann dissipators. We would like to diagonalize the q block L (q) = cL (c) + νL (ν X ) + (cη)L (η X ) = c sin(q/2) D ⊥ − D † ⊥ − (ν/2)r 2 + (cη) cos (q/2)r 2 Each such block produces eigenvalues L (q) |s = −λ q,s |s , that are distinguished by the index s. We are interested in the slowest mode λ q,0 . The NESS is the eigenvector that corresponds to the zero eigenvalue. It belongs to the q=0 block, which results from probability conservation. In the Bloch representation, probability conservation means that 0| L (0) = 0. To obtain the eigenvalues to order q 2 it is enough to Taylor expand the operator to that order. Accordingly, The first term is the zero order term. Here (for X-coupling) it is diagonal in r. For the other coupling schemes it is not necessarily diagonal in r, but for any of them |r = 0 is an eigenstate of the zero-order term.
To find the eigenvalue λ q,0 via perturbation theory one has to sum over different paths that begin and end in r=0. In the case of Eq.(S-38) these paths are composed of hops between near neighbor sites. Second order contributions involve terms with 0| L (q) |r r|L (q) |0 , with r =0. Each transition involves a factor cq or (cη), or (cηq 2 ). Hence only the sites |r| < 2 contribute to the perturbed eigenvalue up to order η 2 q 2 . Furthermore, the (cηq 2 ) transitions are always multiplied by other O(q) transitions, and therefore can be ignored in any second order expansion.
From the above it should be clear that for X-coupling the matrix that should be diagonalized is A convenient way to obtain analytical result is to write the characteristic equation det[λ + L (q) ] = 0 with the above (truncated) matrix, and to substitute an expansion λ q,0 = n a n q n . Then we solve for the coefficients a n iteratively. The outcome is expanded in η to order η 2 . Note that to go beyond second order in η does not makes sense, because the Ohmic master equation and the associated NESS are valid only up to this order.

====== [9] Calculations of A and C
The determination of the A and C coefficients in the main text is obtained analytically from the perturbation theory results of the former section. Numerically we use Eq.(12) with the correlation function of Eq. (20). For the semi-classical B-model, one has the same kernel as in Eq. (16), multiplied by 4 cos 2 (P ). In the latter case, in spite of having zero transition rate for P = π/2, the result for D is finite.