Quantum ergodicity breaking in semi-classical electron transfer dynamics

Can the statistical properties of single-electron transfer events be correctly predicted within a common equilibrium ensemble description? This fundamental in nanoworld question of ergodic behavior is scrutinized within a very basic semi-classical curve-crossing problem. It is shown that in the limit of non-adiabatic electron transfer (weak tunneling) well-described by the Marcus–Levich–Dogonadze (MLD) rate the answer is yes. However, in the limit of the so-called solvent-controlled adiabatic electron transfer, a profound breaking of ergodicity occurs. Namely, a common description based on the ensemble reduced density matrix with an initial equilibrium distribution of the reaction coordinate is not able to reproduce the statistics of single-trajectory events in this seemingly classical regime. For sufficiently large activation barriers, the ensemble survival probability in a state remains nearly exponential with the inverse rate given by the sum of the adiabatic curve crossing (Kramers) time and the inverse MLD rate. In contrast, near to the adiabatic regime, the single-electron survival probability is clearly non-exponential, even though it possesses an exponential tail which agrees well with the ensemble description. Initially, it is well described by a Mittag-Leffler distribution with a fractional rate. Paradoxically, the mean transfer time in this classical on the ensemble level regime is well described by the inverse of the nonadiabatic quantum tunneling rate on a single particle level. An analytical theory is developed which perfectly agrees with stochastic simulations and explains our findings.


Introduction
Statistical physics and physical chemistry normally deal with large ensembles of physical entities. 1 With the advance of single-molecular research on the nanoscale the statistical properties of single separated objects become even more important. In macroscopic systems, fluctuations play an additional subordinated role, away from phase transitions, on top of deterministic mean-field dynamics. In contrast, the dynamics of single elements is fundamentally stochastic, always. Consider a simple two-level system mimicking e.g. electron transfer (ET) between two sites of localization in a molecule, 2-5 as shown in Fig. 1. Within the picture of two diabatic spatially localized electronic quantum states |1i and |2i with energies E 1 and E 2 , an electron can make transitions between these two states due to the tunnel coupling V tun , which provides corrections beyond the adiabatic Born-Oppenheimer approximation (BOP). [2][3][4][5] The electron energies E 1 (x) and E 2 (x) depend on a collective nuclear reaction coordinate x and provide potential curves for the dynamics of this latter one within BOP. Within such a one-dimensional reaction coordinate picture, the rest of the molecular nuclear degrees of freedom and environment act as friction and noise in the dynamics of the reaction coordinate. [2][3][4][5] These friction and noise are related by the fluctuation-dissipation theorem (FDT) at thermal equilibrium. 2,3 Fig. 1 Curve crossing problem in the case of two equal potential curvatures k (i.e. no nuclear frequency change occurs at electronic transitions). Diabatic electron energy levels E 1,2 (x) provide harmonic potentials for the dynamics of nuclear or molecular reaction coordinate x. x 0 is the nuclear equilibrium shift for different electronic terms, and e 0 is the corresponding electron energy difference. l = kx 0 2 /2 is nuclear (molecular) reorganization energy.
The potential minima correspond to the equilibrium nuclei positions in the corresponding electronic states. Note that the spatial localizations of the electronic wave functions c 1,2 (r) = hr|1,2i have nothing in common with the nuclear x. 2,3 The electronic transitions 1 -2 or 2 -1 are accompanied by a spatial jump of the electron, which occurs instantly. 2 This is the so-called curve crossing problem, which provides a basic model for the quantum transport in condensed matter. 2,3,5 Importantly, it applies not only to electron transfer in various molecular systems, 2,6,7 but also to electronic transitions in colloidal semiconductor quantum dots. [8][9][10][11][12][13] Within this model at a single-molecular level, the electron stochastically fluctuates between two sites of localization, say r 1 and r 2 , exhibiting a two-state or dichotomous jumping process in space r. The probability density of the electron is given by |c i (r)| 2 , depending on its quantum state ''i''. On the ensemble level, a standard approach is to derive a kinetic equation for the reduced density matrix r ij (x,v,t) of the ''electron + reaction coordinate'' (v = : x), and one for the electronic subsystem only, P ij ðtÞ ¼ Ð r ij ðx; v; tÞdxdv. [2][3][4][5] In the simplest, ultimately reduced case the electronic populations p i (t) = P ii (t) just follow a Pauli master or balance equation, which corresponds to a classical Markovian two-state process of electron jumps with some quantum rates. Such and similar Markovian kinetic equations describe quantum relaxation dynamics of an ensemble of many identical systems. 2,3 A very profound question is: can such kinetic equations describe the statistics of single-trajectory events? Contrary to a popular belief, which currently dominates the literature, the general answer is surprisingly NO, which is the focus of this paper. It affects a large body of current research on nanophysics and nanochemistry. It will be shown below that within a seemingly classical (on the ensemble level) solvent-controlled ET regime the genuine statistics of single-electron transfer events is very different from the one given by the master equation for electronic populations. This is because ET retains its quantum nature and is intrinsically non-Markovian. The conflict between the ensemble and trajectory descriptions means that ergodicity is fundamentally broken in this regime. These results are very important e.g. for the research on blinking statistics in single colloidal quantum dots [11][12][13] because the discussed kinetics equations are central to the theory approaches based on continuous spectral diffusion. [8][9][10] Basically, our main conclusion is that the standard kinetics equations may not be suitable to describe the statistics of single trajectories in a solvent-or diffusion-controlled regime, even though they can nicely describe the ensemble kinetics at the same time. Here we reveal a profound ergodicity breaking, a deep conflict between single-trajectory and ensemble level descriptions.
The important observation that trajectory and density descriptions of stochastic and kinetic processes can be very different has been made in ref. 14. Subsequently, it has been realized 15-22 that this is a general feature of continuous time random walk (CTRW) dynamics featured by divergent mean residence times (MRTs) in the traps, [23][24][25][26] and the related descriptions such as the fractional Fokker-Planck equation (FFPE). 26 Such FFPEs may fail to describe the statistical properties of single trajectories. 20 For example, subdiffusive FFPEs can be used to derive a general expression for fractional velocity (subvelocity) of biased subdiffusion in tilted periodic potentials. [27][28][29] However, the FFPE is useless to find the mean subvelocity of a single particle, no matter how long its trajectory is. The latter quantity remains always random and exhibits some universal fluctuations. 16,20,29 Averaging over these fluctuations gives the ensemble value of fractional velocity, which can be found from the FFPE. 20 Ergodicity can be understood in various senses. The primary definition is that the time and ensemble averages of a stationary physical observable, say x(t), coincide in the limit of infinitely long trajectories and infinitely large ensembles, correspondingly. 30 This is ergodicity in mean. 30 Clearly, for any finite system both averages are different. They are also trivially different for any nonstationary dynamics. 30 Clearly, if the phase space of a dynamical system is separated into some unconnected domains, the corresponding dynamics cannot be ergodic in principle if the statistical averaging is done on the whole phase space. The reason for ergodicity breaking of CTRW dynamics with divergent MRTs is very different. It is rooted in the absence of stationarity of the increments 22 dx(t|t 0 ) = x(t + t 0 ) À x(t 0 ) with respect to the sliding time point t 0 , which is used for computing a timeaverage. For such processes, the increments are never stationary, even in the limit t 0 -N. Hence, the time average of [dx(t|t 0 )] 2 with respect to t 0 can never coincide with the ensemble average h[dx(t|t 0 )] 2 i, which depends on t 0 because of the absence of stationarity, for any finite t 0 . † This was named as weak ergodicity breaking. 16,31 The relation of weak ergodicity breaking to aging was first realized by Bouchaud. 31 The corresponding systems age for infinite time and never reach a stationary limit. In the real world, this can be of course mostly a transient phenomenon. However, even some standard normal-diffusion models like diffusion in Gaussian disordered potentials with short-range correlations 32,33 can be mesoscopically anomalous [34][35][36] and non-ergodic, 36 even if they are truly ergodic in the limit t 0 -N. In the current literature, the ergodicity of diffusion processes is mostly understood namely in the sense of mean-squared increments. 37 Nevertheless, ergodicity can also be understood in the sense of higher moments, autocorrelation functions, and distribution densities. 30 In the latter case, the question is: 30 can we deduce the probability density of an observable from its single infinitely long trajectory? If yes, the process is ergodic in the probability density sense.
The ergodicity violation which we describe in this paper occurs in the sense of a fundamental conflict between the single trajectory and equilibrium ensemble descriptions. Namely, we ask the question: are the statistics of electronic transitions the same when derived from single electron trajectories (as in single-molecular experiments 38 ) and, alternatively, from a common density description with initial equilibrium distribution of the reaction coordinate (as in the majority of pertinent theories developed thus far and macroscopic experiments)? We show that they are very different in the solvent-controlled ET regime, even if the reaction coordinate x(t) dynamics is fully ergodic. This is because ET remains fundamentally quantum-mechanical and non-Markovian even within this seemingly classical Markovian regime. Hence, standard master or kinetics equations may fail completely in describing the statistics of single trajectories in such and similar regimes. This is a real surprise overlooked thus far, except for a paper by Tang and Marcus, 9 which does not address, however, the issue of ergodicity. We develop both an analytical theory and stochastic numerics, which remarkably agree, confirming our major results and conclusions.

Theory
2.1 Landau-Zener tunneling, semi-classical curve crossing problem, and Zusman equations Consider a curve-crossing problem in Fig. 1 for two parabolas, E i (x) = k(x À x 0 d 2,i ) 2 /2 À e 0 d 2,i , of equal curvature k (i.e. no nuclear frequency change occurs at electronic transition, for simplicity).
x 0 is the reaction coordinate shift, e 0 is the equilibrium difference of electron energies, and d 2,i is the Kronecker symbol. The quantum system is characterized by the Hamiltonian Ĥ (x) = E 1 (x)|1ih1| + E 2 (x)|2ih2| + V tun (|1ih2| + |2ih1|), and the reaction coordinate x(t) is treated classically. Tunnel coupling V tun is assumed to be coordinate independent, as in most ET theories, which is known as Condon approximation. 2,3 This, of course, does not exhaust all the possibilities. Non-Condon effects in general and, in particular, for conical intersections, where V tun p x À x*, have attracted increasing attention. 39,40 They are beyond the scope of this work, where we focus on the simplest possible minimal model. What is the probability of electronic transition between two localized quantum states |1i and |2i depending on the reaction coordinate velocity v? Within the Landau-Zener-Stückelberg theory (LZS) [41][42][43] applied to this problem the answer is 2 This result presents a milestone in the theory of quantum transport. Here, is the lowest second order approximation in the tunnel coupling V tun , and DE(x) = E 1 (x) À E 2 (x) is the difference of electron energies. The latter result follows immediately from Fermi's Golden Rule quantum transition rate applied at the level crossing point x*, DE(x*) = 0. Here, d(x) is Dirac's delta-function. In the present case, DE(x) = e 0 À l + 2lx/x 0 , where l = kx 0 2 /2 is the nuclear reorganization energy. Once again, the electron tunnel distance has nothing in common with x 0 . Electrons tunnel in space once a transition |1i -|2i, or |2i -|1i takes place. Very importantly: while electronic transitions always show a discontinuous jump-like quantum character, the reaction coordinate dynamics remain continuous. This is the ultimate reason for the quantum breaking of ergodicity, see below. Likewise, blinking of a quantum dot occurs from a light emitting quantum state upon a radiative quantum transition. 8 Depending on the coupling strength V tun and the velocity v = : x at the crossing point P LZ (v) can vary from P LZ (v) E f (v) p |V tun | 2 /|v| (nonadiabatic transition) to one (adiabatic transition).
Within a classical treatment of the reaction coordinate x, one considers it as a particle of mass M subjected to viscous frictional force Zv, with a friction coefficient Z, and zero-mean white Gaussian thermal noise of the environment x(t) at temperature T. Friction and noise are related by the fluctuationdissipation relation hx(t)x(t 0 )i = 2k B TZd(t À t 0 ), where hÁ Á Ái denotes the ensemble averaging. Stochastic dynamics of x(t) follows the Langevin equation: which depends on the quantum state |ii. The electron-reaction coordinate dynamics can be described in a semi-classical approximation by mixed quantum-classical dynamics of the reduced density matrix r ij (x,v,t), where the quantum degree follows quantum dynamics while the dynamics in (x,v) phase space for a fixed quantum state i is classical. Generally, the phase-space dynamics is described by the Kramers-Fokker-Planck equation (KFPE). In the overdamped case, Z ) ffiffiffiffiffiffiffiffi Mk p , the reaction coordinate velocity is thermally distributed, The corresponding semi-classical description is well known under the label of Zusman equations. 4,44 Within it, the reduced dynamics of population densities p i ðx; tÞ : where G(x) is the Golden Rule expression in eqn (3). This equation is obtained from a full reduced density matrix equation description after excluding (projecting out) the dynamics of quantum coherences, and within the so-called contact approximation. 44 One must stress that no rigorous second order approximation in the tunnel coupling V tun is done in eqn (5), but rather nonlocality and memory effects are neglected. 45 The appearance of G(x) p |V tun | 2 in eqn (5) is due to a singular limit of overdamped dynamics. Indeed, since v T -N in this limit, one can safely approximate P LZ (v) E f (v), and also the contact approximation (locality in x) is well justified from the quantum uncertainty relation. 4 We confirm below this very important remark while doing stochastic simulations of electron trajectories. In these simulations, we use generically P LZ (v), instead of f (v), which, however, does not change the results for overdamped dynamics. For a strong electron-nuclear coupling (l c V tun ), in the absence of inertial effects, and in the limit where the quantum effects in the reaction coordinate dynamics are entirely neglected, this contact approximation is well justified. 4,44 It presents a very important reference point, which allows also for further generalizations toward anomalous subdiffusive dynamics of the reaction coordinate. 10 Indeed, within this approximation one obtains very elegant and important analytical results. Consider first very small V tun , with the reaction coordinate being thermally equilibrated, p is a thermal width. Then, the nonadiabatic quantum transition rate is with activation energies E i on e 0 is famously known as a Marcus parabola. Note that in this respect the so-called inverted regime of electron transfer for e 0 4 l is entirely a quantum-mechanical feature which is physically impossible within an adiabatic classical treatment.
With the increase of V tun the reaction coordinate dynamics become even more important and can limit the overall rate. Assuming that the reaction coordinate is thermally equilibrated in the initial electronic state, the following master equation : has been derived 4,44,45,50-52 from eqn (5) with the rates In eqn (8) is the mean escape time in the parabolic potential with cusp, 45 and t = Z/k is the reaction coordinate relaxation time, or the Debye solvent correlation time. Here, 2 F 2 (a,b;c,d;z) is a generalized hypergeometric series. 53 It is worth noting that this theory is restricted by the parameter domain V tun { k B T, l. The relaxation of populations is single-exponential with p 1,2 (N) = 1/[1 + exp(AEe 0 /k B T)], and k = k 1 + k 2 . Physically, this requires sufficiently large activation barriers. For 44 Hence, for such large activation barriers, eqn (8) takes the form which is known as the Zusman rate, 44 where is an adiabaticity factor. Moreover, for k ad c 1 and e 0 o l which is the adiabatic Marcus rate. For a particular case 1,2 coincides with the Kramers rate for the adiabatic transitions in the cusp potential 33 consisting of two pieces of diabatic curves in Fig. 1, and Hence, for a sufficiently large V tun or sufficiently large t (solvent control), so that k ad c 1, ET becomes purely classical and adiabatic within this ensemble description. This is the so-called solvent-controlled adiabatic ET. The rate expressions (6), (11), (12) and (14) are traditionally used to interpret and analyze the experimental data. 7 The fact that a change of the solvent relaxation time t can control the rate of intramolecular ET has been shown e.g. in ref. 54. Long-range electron transfer occurs typically in the non-adiabatic regime, while intramolecular electron self-exchange processes in molecular compounds and mixed-valence systems can occur adiabatically. 1,2,6 3 Results and discussion

Trajectory description
In this work we focus our attention on a stochastic trajectory counterpart of this well-known ensemble theory. It can be obtained as follows, quite in the spirit of a surface hopping approach. 39,55 We propagate overdamped (with M = 0) Langevin dynamics (4) on one potential surface. Once the threshold x* is reached quantum hopping onto another surface occurs with the LZS probability (1), where v = dx/dt, dt is the time integration step, and dx is the x displacement by crossing the threshold. After a quantum jump, Langevin dynamics are continuously propagated on another surface, until the next jump occurs. Note that even if for dt -0 the formal limit of dx/dt does not exist in a mean-square sense for the strictly overdamped dynamics, at any finite dt, v is finite. The overdamped dynamics of the reaction coordinate lead, however, to effective linearization of eqn (1) the results do not depend on whether we use eqn (1), or (2) in simulations of the overdamped dynamics. This is our first remarkable result which is completely confirmed by numerics and agrees with the Zusman equation theory. We consider the symmetric case e 0 = 0 further in this work. Numerical simulations of eqn (4) in a nondimensionalized form ‡ for a fixed electronic state were performed using the stochastic Heun method 56 with a fixed time integration step dt, which has been varied to achieve the convergency of numerical results. In most simulations, dt = 10 À5 . For the strictly overdamped (M -0) dynamics and P LZ (v) E f (v), the dynamics of an ensemble of particles based on single-trajectory simulations corresponds precisely to the density dynamics in eqn (7). By propagating many (10 4 ) particles simultaneously starting from the quantum state ''1'' and distributing initial x(0) in accordance with P (eq) 1 (x), we can keep track of the state populations. The corresponding results in Fig. 2 for a sufficiently high activation barrier l/4 of 2.5k B T, which is typical for many molecular systems, e.g., for ET in the azurin dimer, 57 agree remarkably well with the theoretical result in eqn (6)- (10). In other words, the result of ensemble-averaged trajectories nicely agrees with the analytical solution of Zusman equations. The ensemble kinetics is practically single-exponential. For a very small V tun , ET is non-adiabatic and characterized by the MLD rate. Upon the increase of V tun , the adiabatic transport regime is gradually established. It is almost achieved for V tun = 0.03, as shown in Fig. 2.

Statistics of single trajectories from a master equation perspective
Single-trajectory electron transition statistics is commonly characterized by residence time distributions (RTDs) in the electronic states. Such a RTD at the ensemble level can be obtained by preparing all the particles in one state, with the reaction coordinate initially thermally equilibrated and taking out particles once they jump to another state, until no particles remain in the initial state. This is how one standardly proceeds within a master equation theory. For the considered high activation barrier, the corresponding survival probability F 1 (t) decays single-exponentially, see Fig. 3, however, with the rate G 1 , which is different from the above k 1 . Indeed, by repeating straightforwardly an analytical derivation in ref. 45 in place of eqn (8). Another derivation of this result is given in Appendix A. This result has a very simple interpretation. Namely, the average time to make a transition is the sum of the average time to reach the threshold x* and the inverse of the nonadiabatic tunneling rate. Indeed, numerics remarkably agree with this simple and insightful result based on the density description, see Fig. 3.

Non
-Markovian yet exponential kinetics. Trajectory simulations contain, however, in principle much more information than a density description can deliver. In particular, we can study the residence time distributions (RTDs) in the electronic states directly by observing a very long single trajectory. We proceed further by noticing that for Markovian dynamics it must be G 1,2 = k 1,2 . This is indeed the case in the nonadiabatic ET regime. However, the dynamics of single-electron transitions become increasingly non-Markovian upon taking adiabatic corrections into account with the increase of V tun . This is in spite of a single-exponential character of the ET kinetics at the ensemble level! Ref. 58 already pointed out that in a similar very paradoxical situation: a highly non-Markovian bursting process can have a nearly exponentially decaying autocorrelation function. Indeed, for most researchers the exponential decay of survival probability   Fig. 3 would imply that electronic state populations undergo a two-state Markovian stochastic process with two equal rates G 1 = G 2 . If so, then on average there would occur only one electronic transition within a mean residence time interval 1/G 1 . However, a short inspection of a single trajectory realization of electronic transitions in a near to adiabatic regime depicted in Fig. 4 reveals immediately that this is not the case. Many transitions can occur on the time scale 1/G 1 . This reveals clearly a non-Markovian character of electronic transitions, which remains hidden in the traditional ensemble theories operating using Markovian master equations for the electronic populations, like eqn (7). Bursting provides such a visual proof. 58 Hence, the actual statistics of electronic transitions are expected to be very different from those implied by the ensemble kinetics measurements. This implies that ergodicity is broken. Note that a popular statement that in an adiabatic ET regime an electron just follows nuclear rearrangement is in fact rather misleading on the level of single electron trajectories. This is so because electrons jump immediately at the level crossing (in the contact approximation) and not after nuclei complete their rearrangement. ET remains quantum even within this adiabatic seemingly fully classical regime! And namely this causes the quantum breaking of ergodicity described next.

Actual statistics of single trajectories
Indeed, the study of survival probabilities based on very long single trajectories surprisingly indicates that ergodicity is broken in this profoundly non-Markovian regime. The corresponding survival probability in a state is depicted in Fig. 5(a). It is profoundly nonexponential, very different from the corresponding ensemble result in Fig. 3. The rate G 1 describes only the tail of distribution, which is the initially stretched exponential. It can possess also an intermediate power law regime for a larger V tun , see part (b) in Fig. 5, where the exponential tail has a weight of less than 10%. Very surprisingly, the mean residence time is well described by the inverse of the Marcus-Levich-Dogonadze rate, ht i i = 1/k (nad) i , at the level of single trajectories. This result must be contrasted with the prediction of the equilibrium ensemble theory, ht i i ens = 1/G i , extended ad hoc onto single trajectories. The latter one is in fact completely wrong because the electron jump process is profoundly non-Markovian. Indeed, to derive the correct distribution of singleelectron residence times from Zusman eqn (5) one must solve a very different initial-value boundary problem. 8 Indeed, after each jump the electron starts its evolution at a very non-equilibrium value of the reaction coordinate x E x*. For example, in state ''1'', p 1 (x,0) = d(x À x* + e), with e -0, and a radiative boundary condition is placed exactly at x = x*. An analytical solution of this problem in the Laplace space is given in Appendix A. The corresponding theoretical result in eqn (A18), which is different from the one attempted by Tang and Marcus,8 nicely agrees with the numerics in Fig. 5. Moreover, a simple and nice analytical approximation follows from eqn (A13), in the adiabatic limit k (nad) i t (ad) i c 1: where E 1=2 À ffiffi ffi z p ð Þ¼e z erfc ffiffi ffi z p ð Þ is the Mittag-Leffler function of the fractional order 1/2 and argument À ffiffi ffi z p , and  is a fractional transfer rate, see Appendix A for details. For high activation energies, E (a) i c k B T, This approximation excellently describes the non-exponential transport regime in Fig. 5, with a i E 1, which covers about 90% of the probability transfer in Fig. 5(b). In accordance with it, the initial decay of the survival probability is always a stretchedexponential, Here we introduced a characteristic time The corresponding residence time density is Note that even though c i (t) does contain a 1 ffiffi t p scaling part, our result is very different from the one by Tang and Marcus. 8,9 First, it describes a stretched-exponential and not a power law distribution. Second, our t i,+ in (19) is also different from the corresponding critical time in ref. 8  and c i (t) p t À3/2 , which ends with an exponential tail, which is described by the rate G i in eqn (15). This power law regime and an exponential tail were obtained in ref. 8. However, the exponential tails here and in ref. 8 and 9 have different rates. All in all, the corresponding statistics are generally very different from the one following from the custom equilibrium ensemble description of adiabatic electron transfer. Strikingly enough, while the ensemble kinetics is strictly exponential, the single-electron transfer statistics is described by the Mittag-Leffler relaxation kinetics in the strictly adiabatic ET regime. Therefore, one can speak about profound ergodicity breaking in this context. The observed striking result, namely that the mean residence time is always given by the inverse MLD rate, is explained within our analytical expression (A9) for the survival probability derived in Appendix A. Moreover, our analytical result, which corrects the earlier one by Tang and Marcus,8 explains also why the tail of distribution is always given by the Zusman type rate G i in eqn (15), which becomes the Marcus adiabatic rate in the strict adiabatic limit. Most of these remarkable features are missed in the previous analysis by Tang and Marcus. 8,9 Moreover, the result that the mean residence time is always the inverse of a non-adiabatic MLD rate can also be explained within a modification of the classical level-crossing theory. 30 For this, let us take formally into account small inertial effects (keeping M finite). Then, the velocity process v(t) is not singular anymore. Consider dynamics in the state i. Assuming stationarity of x(t), the averaged number of level crossings n i (T) within a very long time interval T is according to ref. 30 n i (T) = TP (eq) i (x*) h|v(t)|i x(t)=x* , and hence, t i h i À1 ¼ lim the same token and taking into account probability (1) to make a quantum jump to another state at each level crossing, we obtain Indeed, consider an auxiliary variable j(t) = x(t) À x*. The crossings occur at random times t i , and dðjðtÞÞ ¼ P Note that if the transitions occur at each crossing, the integral of l.h.s. of eqn (22) within a time-interval T would give the number of transitions to another electronic surface within this time interval. Since the probability of such transitions is merely P LZ (v(t i )), one should consider instead of eqn (22). Then, the time-average of z(t) in the case of infinitely long T yields the inverse mean time between the switching events, or the mean transition rate. Assuming ergodicity of the point process t i , the time-average can be replaced with the ensemble average using P(x,v) = P (eq) (x)P M (v), which yields (21). Note a very paradoxical situation: even if z(t) is an ergodic process, the process of electron fluctuations between two localization sites is generally nonergodic, whenever its equilibrium ensemble and single trajectory statistics are in conflict. Averaging in (21) with Maxwellian equilibrium P M (v) yields an important result where is a renormalization function taking inertial effects into account. It is expressed via a Meijer G-function, 53 and v 0 = p|V tun | 2 x 0 /( hl) is a characteristic tunnel velocity. Numerically, R(z) E exp(À1.57z 0.9 ) for 0 o z o 0.1 with the accuracy of about 10%. In the formal overdamped limit, lim , and we obtain ht 1 i À1 = k (nad) 1 , in agreement with numerics. Moreover, we also did numerics which include the inertial effects in eqn (4), within the contact approximation, and confirm the analytical result in (24) and (25), see Fig. 6. The numerical results for the renormalization function R(z) (symbols) are compared therein with the analytical result in eqn (25). The agreement is remarkable indeed, and a detailed treatment of the inertial case will be presented elsewhere. The observed ergodicity breaking is thus not an artifact of the overdamped singular approximation. It expresses the profoundly quantum nature of electron transfer even in the adiabatic regime as manifested at the level of single-trajectory dynamics.

Conclusions and outlook
As a major result of this work, equations like Zusman equations and other semi-classical and even fully quantum-mechanical ensemble descriptions should be used with great care to describe the properties of single electron trajectories in a profoundly non-Markovian ET regime, which can mistakenly be identified as a Markovian one (due to single-exponential kinetics) on the ensemble level. The kinetics equations involving the reaction coordinated dynamics must be used for this with a very nonequilibrium initial preparation of the reaction coordinate, which reflects the quantum nature of electron transfer, always, even within a seemingly purely classical adiabatic regime. Tang and Marcus earlier observed this striking feature. 9 However, their attention was not drawn to it as a manifestation of profound nonergodic effects, which are currently in the limelight, and they missed several very important features we discussed.
Single-trajectory events can display a stretched exponential and power law statistics in the case of strictly exponential ensemble kinetics. This effect can be relevant e.g. for adiabatic ET in mixed-valence systems and molecular compounds. 1,6 Even long-range ET in proteins, with V tun exponentially decaying with the tunneling distance, can enter the corresponding transfer regime when the reaction coordinate relaxation time is in the nano-and millisecond range as e.g. in cytochrome c, where t B 200 ns. 7 Another class of important experimental systems is provided by colloidal quantum dots such as CdSe-CdS nanocrystals, 12,13 and quantum dot-molecule complexes. 59 It is very important that already a simplest ET model of textbook character and standard use 1-3 studied in this work reveals the effect. Of course, many important issues remain to be addressed and clarified in the future. These are, in particular: (1) will the effect persist also for crossing multi-dimensional electron energy surfaces? (2) What is the role of non-Condon effects, e.g. for conical intersections? (3) What is the role of long-range memory effects in the reaction coordinate dynamics? The latter ones were revealed e.g. experimentally for ET in protein systems. 38 The answer for question (1) seems to be intuitively clear: YES, the effect should persist whenever diffusion of the reaction coordinate on the electron energy surface limits the overall ET rate within the ensemble description. However, additional research is nevertheless welcome. Many more questions can, of course, be raised and addressed in follow-up research.
Beyond doubt, such a nonergodic behavior is pertinent also for diffusion-controlled dynamics in the case of low activation barriers, where power law features emerge already on the ensemble level. This can also be relevant e.g. for the blinking statistics in quantum dots in non-exponential regimes, whenever the reaction coordinate dynamics is very essential. 10 In this case, the predictions based on equilibrium density dynamics are generally inappropriate to describe the actual statistics of single trajectory events. This is especially true for anomalously slow subdiffusive dynamics with memory, which is the subject of a separate follow-up work. In the case of normal reaction coordinate diffusion, a non-equilibrium ensemble perspective does deliver the correct result for single trajectories, as we showed in this work. Will it be also true for anomalously slow stochastic dynamics with long-range memory? There exist profound doubts that it is so. The intrigue remains because the formal solution of the corresponding non-Markovian Fokker-Planck equation of viscoelastic subdiffusion 10,60 can be very different from a more physical treatment of the underlying stochastic dynamics with memory, which properly takes the slow modes of the environment into account, as explained in detail in ref. 60. This problem will be addressed in a followup work. The discovered non-ergodic features in a simple and well-known model of normal-diffusion charge transport dynamics are expected to influence a large body of current research on nanoscale transport phenomena. The author is confident that many other ensemble transport models should be reexamined from a trajectory point of view with focus on possible non-ergodic features.

A Appendix
Let us consider first electrons starting in an electronic state ''i'' at t = 0 with the initial equilibrium nuclei distribution p (eq) (x). The subindex ''i'' is dropped for a while for simplicity of notations. Given a singular sink term S(x) = v 0 d(x À x*), with v 0 = p|V tun | 2 x 0 /( hl), the probability density p(x,t) develops in time as follows: 61 pðx; tÞ ¼ p ðeqÞ ðxÞ À v 0 ð t 0 dt 0 G x; t À t 0 jx Ã ð Þ p x Ã ; t 0 ð Þ; (A1) Fig. 6 Renormalization factor R as a function of is angular frequency of the reaction coordinate. The theoretical result (full line) in eqn (25) is compared with the numerical results (symbols) obtained for V tun = 0.04, l = 800, and T = 0.1, and varying z through o 0 .
From it,p Next, by integrating eqn (A3) using the normalization of probability, we obtaiñ The first line in eqn (A5) says that the RTD is c(t) = v 0 p(x*,t) because, quite generally, F(s) = [1 Àc(s)]/s. Next, we note that in the limit s -0, i.e. G(x*,s|x*) has a singularity at s = 0, which is quite obvious because lim t!1 G x Ã ; t x Ã j ð Þ¼p ðeqÞ x Ã ð Þ. Here lies the major difference of our treatment with the one in ref. 8, which we correct thereby. Subsequently, we use the decomposition G(x*,s|x*) À p (eq) (x*)/s = t(s), wherẽ tðsÞ ¼ and note thatt(0) = t (ad) , see ref. 45, which yields eqn (9). Furthermore, because v 0 p (eq) (x*) = k (nad) , we finally havẽ (A8) where we restored the subindex of the electronic state. Since F i (0) = ht i i is the mean residence time in the state ''i'' and G i = 1/ht i i is the corresponding rate in the exponential regime of high barriers, we obtain immediately from eqn (A8) our result in eqn (15). Next, let us consider the particles starting at some x 0 for t = 0. Then, in eqn (A1) one must replace p (eq) (x) with G(x,t|x 0 ). Repeating the above derivation for x 0 = x*, we obtaiñ From it, two of our central results follow immediately. First, in this case ht i i = 1/k (nad)