{'fraction_non_alphanumeric': 0.037260734315751234, 'fraction_numerical': 0.011375763493108693, 'mean_word_length': 4.38476586379753, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dissipative Two-Electron Transfer

We investigate non-equilibrium two-electron transfer in a model redox system represented by a two-site extended Hubbard model and embedded in a dissipative environment. The influence of the electron-electron interactions and the coupling to a dissipative bosonic bath on the electron transfer is studied in different temperature regimes. At high temperatures Marcus transfer rates are evaluated and at low temperatures, we calculate equilibrium and non-equilibrium population probabilities of the donor and acceptor with the non-perturbative Numerical Renormalization Group approach. We obtain the non-equilibrium dynamics of the system prepared in an initial state of two electrons at the donor site and identify conditions under which the electron transfer involves one concerted two-electron step or two sequential single-electron steps. The rates of the sequential transfer depend non-monotonically on the difference between the inter-site and on-site Coulomb interaction which become renormalized in the presence of the bosonic bath. If this difference is much larger than the hopping matrix element, the temperature as well as the reorganization energy, simultaneous transfer of both electrons between donor and acceptor can be observed.


I. INTRODUCTION
Electron transfer is a key process in chemistry, physics and biology 1,2,3,4 encountered in, e.g., chemical redox processes, charge transfer in semiconductors and the primary steps of photosynthesis. In condensed polar environments the process involves strong coupling to the underlying nuclear motion and is usually dominated by the nuclear reorganization that accompanies the charge rearrangement. A quantum mechanical description of electron transfer in such a dissipative environment is given by the spin-boson model 5,6 and its variants; this model accounts for the essential energetics and dynamics of the process, such as the non-monotonic dependence of the transfer rate on the energy asymmetry, the energy difference between the initial and final electronic states.
Although standard descriptions of such processes focus on single-electron transfer 1,4,5,6 , two-electron transfer has been suggested as the dominant mechanism in some bioenergetic processes that occur in proteins 7,8 , transfer in transition metal complexes 9,10 , electrode reactions 11 , artificial photosynthesis and photoinduced energy-and electron-transfer processes 12 , biological electron transfer chains 13 , transfer in fuel cells 14 as well as in DNA 15 . Further examples are selfexchange reactions such as Tl(I)/Tl(III) and Pt(II)/Pt(IV) 16 and electronpair tunneling 17,18,19 in molecular electronic devices.
The theoretical description of two-electron transfer dynamics differs fundamentally from its singleelectron counterpart. More than two states have to be considered 20,21 and electron correlations induced by the Coulomb repulsion and the coupling to the environment need to be accounted for. Usually, the on-site Coulomb interaction in molecules is much larger than the intersite interaction. 22,23,24 However, due to the polarization of the local environment the short range interaction may be strongly screened. Then, the inter-site interaction V can be of the same order or even exceed the on-site Coulomb interaction U . 24,25 While U favors a homogeneous charge distribution, the inter-site interaction V inclines spatially inhomogeneous charge accumulation. Since the non-equilibrium dynamics is governed by the energy difference U − V , the competition between both interactions influences strongly the type of charge transfer dynamics. Depending on the sign of the energy difference a single concerted two-electron step or two sequential single-electron steps may occur.
In this paper we consider a system comprised of a donor (D) and an acceptor (A) site. They share two electrons which are coupled to a non-interacting bosonic bath. Such a donor-acceptor system has four different states: two doubly occupied donor (D 2− A) and acceptor (DA 2− ) states, and two degenerate states D − A − with one electron each on the donor and acceptor site (with different spin). Their energy difference depends on the difference between on-site and inter-site Coulomb repulsion as well as the bias ε which we do not consider here. The transition D 2− A → DA 2− occurs as a concerted transfer of two electrons or an uncorrelated sequence of one-electron transfer events during which the intermediate D − A − is formed. The transfer rate of each electron may be different and shows a non-monotonic behavior on the energy asymmetry between the states. In this paper, we are mapping conditions under which the system performs concerted two-electron transfer or a sequential single-electron process. To this end we study the non-equilibrium dynamics of the donor-acceptor system initially prepared with two electrons at the donor site. We evaluate the rates for single-electron transitions and an electron-pair transfer in different regimes of the Coulomb repulsion and environmental response.
The occurrence of such a correlated electron-pair transfer can be already understood within a donor-acceptor system, decoupled from the environment in which the strong on-site Coulomb repulsion exceeding considerably the inter-site repulsion. We start from a doubly occupied excited donor state D 2− A compared to the D − A − ground states. Energy conservation implies a concerted electron transfer. If the transfer matrix element ∆ is much smaller than this energy difference the transfer occurs as a tunneling process of an electron pair in which the intermediate states D − A − are occupied only virtually analogous to a "superexchange" process (see, e.g., Ref. 1).
In the present paper we investigate the effect of coupling to a dissipative bosonic environment, with a total number of two electrons occupying donor and acceptor sites. These two electrons experience the onsite Coulomb repulsion U when occupying the same site and the Coulomb repulsion V when occupying different sites. In this paper we restrict ourselves to the simplest case where donor and acceptor are each modeled by a single molecular orbital. In such a system the dif-ferenceŨ = U − V is crucial for the dynamics. The coupling to the bosonic bath has two major effects: (i) the renormalization 26,27 of the on-site Coulomb repul-sionŨ toŨ eff and (ii) dephasing as well as dissipation of the energy from the donor-acceptor system to the bath. The latter leads to the damping of coherent oscillations that would otherwise exist between the quantum states of the related molecule and, beyond a characteristic coupling strength, to incoherent dynamics of the electron transfer process. These considerations lead us to a dissipative two-site Hubbard model, a minimal model that captures the essential physics comprising correlations between electrons and their coupling to the dissipative environment. It is discussed in detail in Sec. II. For a comparison to experimental results it has to be supplemented by ab initio calculations of the parameters.
The equilibrium properties of the model have been previously studied 28 using the numerical renormalization group (NRG), and the real-time dynamics has been investigated 29,30 using a Monte-Carlo technique at high temperatures where only incoherent transfer is present. In these Monte-Carlo calculations, the effective Coulomb interaction was chosen to beŨ eff > 0 and no electronpair transfer has been reported. Two-electron transfer in a classical bath has been discussed in Ref. [20] in the framework of three parabolic potential surfaces (for the four states D 2− A, D − A − and DA 2− ) as a function of a single reaction coordinate. A generalization to donorbridge acceptor systems is given in Ref. [8,31].
Although two-electron transfer was observed in some regimes of system parameters in the high-temperature limit, considering a classical bath, it seems reasonable to expect that, at least between identical centers, electronpair tunneling processes are particularly important at temperatures corresponding to energies smaller than the effective energy difference between initial and intermediate statesŨ eff . At these temperatures single-electron transfer cannot be activated (see Section VI). Therefore, we focus on the low-temperature regime where the transfer is dominated by nuclear tunneling and where the bosonic bath has to be treated quantum mechanically. Due to the nuclear tunneling the electron transfer rate is constant over a wide temperature range from zero temperature up to temperatures where thermal activation becomes more important 32 . In this low-temperature regime, we employ the time-dependent NRG 33,34,43 (TD-NRG) which covers the whole parameter space from weak to strong dissipation. The NRG is an accurate approach to calculate thermodynamics and dynamical properties of quantum impurity models. 35,36,37,38 For further details of the NRG we refer to the recent review 39 on this method.
The paper is organized as follows. In Sec. II we introduce the model. Its high temperature behavior obtained from the Marcus theory is described in Sec. III. Section IV introduces the NRG method, its extension to non-equilibrium and its application to the present problem. In order to gain a better understanding of the nonequilibrium dynamics presented in Sec. VI, we summarize the equilibrium properties of the model in the Sec. V. We present a detailed discussion of the real-time dynamics in Sec. VI. Therein, we focus on the time evolution of occupation probabilities of the different electronic states as the key observables. In particular, when the dynamics can be described in terms of rate processes, the dependence of the single and electron-pair rate on the Coulomb repulsion parameters is analyzed. A summary of our results is given in Sec. VII.

II. MODEL
We consider a model of a two-electron/two-site system coupled to a bosonic bath. It is defined by the Hamiltonian with where c iσ and c † iσ denote annihilation and creation operators for fermions with spin σ on the donor (i = D) and acceptor (i = A) sites. The Hamiltonian H el corresponds to an extended two-site Hubbard model, with on-site energies ε i , hopping matrix element ∆, on-site Coulomb repulsion U and an inter-site Coulomb repulsion V between one electron on the donor and one electron on the acceptor. The differenceŨ = U − V measures the excess energy needed to move an electron between the two sites. Such a two-site Hubbard model without coupling to a bosonic bath has been investigated in the context of electron transfer in Ref. 40.
The Hamiltonian H b models the free bosonic bath, with boson creation and annihilation operators b † n and b n , respectively. The electron-boson coupling term, H coupl , has the standard polaron form with the coupling constant for donor and acceptor given by g D λ n and g A λ n , respectively. In what follows we set ǫ D = −ǫ A = ε 2 and g A = −g D = 1. The latter choice implies that the polar bath is coupled to the change in the electronic density This two-site electron-boson Hamiltonian conserves the number of electrons iσ c † iσ c iσ and the square of the total spin S 2 as well as its z-component S z . The Hilbert space can therefore be divided into different subspaces. In the subspace with one electron and S z = 1/2, the model is equivalent to the spin-boson model 28 . Here, we consider the subspace with two electrons and S z = 0 which is spanned by the states |1 = | ↑↓, 0 , |2 = | ↓, ↑ , |3 = | ↑ , ↓ , and |4 = |0, ↑↓ with the notation |A, D describing the occupation at the donor (D) and acceptor (A) sites. The four-dimensional basis in the two-electron subspace is displayed in Fig. 1. We define the following observableŝ which measure the doubly occupancyd D (d A ) on the donor (acceptor) site andn DA the combined population of the states | ↑, ↓ and | ↓, ↑ . Note that in some works the states | ↑↓, 0 and |0, ↓↑ are referred to as localized states 29 . We call them doubly occupied states while the term localization is used below for the self-trapping mechanism.
Consider the 4×4 Hamiltonian matrix in the electronic subspace (M ) ij = i|H|j (i, j = 1, . . . , 4). Introducing the notation and shifting the Hamiltonian by a constant V leads to Therefore, the dynamics of the system is governed by the energy differenceŨ which replaces the on-site Coulomb repulsion U . If screening of the local Coulomb repulsion U 24,25 is sufficiently large, U changes its sign and become effectively attractive. A large inter-site Coulomb repulsion V favors an inhomogeneous charge distribution. It is convenient to rewrite the diagonal matrix elements of the doubly occupied states in the form and Compared with the matrix elements of states correspond- we can easily see that the electron-boson coupling generates an effective renormalized interactioñ The renormalized interactionŨ eff determines the energy difference between D 2− A (DA 2− ) and D − A − , and constitutes the only Coulomb interaction parameter in the present model. The renormalization stems from a bosoninduced effective electron-electron interaction, already familiar from the Holstein model 26 . Note that an artificial energy shift is present in the single-electron subspace (spin-boson model) 6 , however the two states | ↑ and | ↓ are shifted in the same direction which can be handled by resetting the zero of energy.
In analogy to the spin-boson model 5,6 , the coupling of the electrons to the bath degrees of freedom is completely specified by the bath spectral function The spectral function characterizes the bath and the system-bath coupling, and can be related to the classical reorganization energy 6 (classical in terms of boson degrees of freedom) which measures the energy relaxation that follows a sudden electronic transition. The oneelectron transfer and the correlated two-electron transfer are associated with reorganization energies E α1 and E α2 , respectively. For a single-electron transfer, e.g., and the corresponding energy for a correlated two- The model we are considering here is completely specified by the parameters ∆, α,Ũ, ε and the bosonic spectral function. In the molecular electron transfer problem the latter function reflects intramolecular vibrations and the solvent (e.g., water or protein) or environment. Its solvent component can be estimated from the solvent dielectric properties or a classical molecular dynamics simulation. In the present paper we assume an Ohmic bath model: with a cut-off at energy ω c . This choice yields the reorganization energy E α1 = 2αω c and the energy shift U eff =Ũ − 2αω c . All parameters and physical quantities are defined in units of ω c . Its order of magnitude for the intermolecular mode spectrum of a polar solvent is 0.1eV.

III. THE HIGH-TEMPERATURE LIMIT: MARCUS THEORY
In the high-temperature, limit electron transfer is usually described using Marcus theory 4 as a rate process within classical transition state theory. Extensions that take into account the quantum nature of the nuclear motion in the weak electronic coupling limit (the so called non-adiabatic limit) are also available 4 , however for simplicity we limit ourselves in what follows to the classical Marcus description. The Marcus electron transfer rate can be evaluated for any amount of transferred charge: the latter just determines the renormalized potential surface parameters that enter the rate expression. Singleelectron transition rates are given by In the case ∆ ≪ |U eff | second-order processes are possible that involve only virtual occupations of the states D − A − leading to rates for an electron pair.
The interplay between sequential and concerted twoelectron transfer (in the limit of a classical bath with a single mode or a single reaction coordinate) can be seen from these expressions. In the following we restrict ourselves to the symmetric case (ε = 0). Starting with the initial state D 2− A, we expect concerted twoelectron transfer in the Marcus regime when the rate k pair [D 2− A→DA 2− ] is larger than the rate k single [D 2− A→D − A − ] of the first step of the sequential process which is the case when |Ũ eff | ≫ T and |Ũ eff | ≫ E α1 as well as E α1 ≤ T .
In a parameter region where sequential transfer dominates the rates k single as well as the corresponding backward rates show a nonmonotonic behavior and an inverted regime dependent on the effective Coulomb interactionŨ eff (see Fig. 2).
For incoherent transfer processes (which may happen at large temperatures and for a strong coupling to the bosonic bath), a description of the population dynamics by kinetic equations determined by the rates is given bẏ where d D and d A are the probabilities to have two electrons on the donor and acceptor, respectively. n DA is the combined population of the states | ↑, ↓ and | ↓, ↑ . For the initial condition d D (t = 0) = 1 we obtain n |↑,↓ = n |↓,↑ . For the unbiased Hamiltonian (ε = 0), These kinetic equations can be solved in the hightemperature regime using the Marcus rates from above. In Sec.VI, we have used these equations to extract the low-temperature transition rates by fitting the nonequilibrium dynamics of d D (t), d A (t) and n DA (t) calculated in the incoherent regime with the time-dependent NRG.
For t → ∞ the equilibrium states d A eq , d D eq and n DA eq are reached, where d D eq = d A eq . It follows that dD eq which is according to the Marcus rates dD eq nDA eq = eŨ eff /T . Therefore, in the classical limit we arrive at With the help of the kinetic equations we can describe concerted two-electron transfer, a purely sequential single-electron transfer as well as a combined process which shows first a pair transfer which is followed by a single-electron transfer. As long as the single-electron transfer rates are small (k single and d D eq = d A eq ≈ 0.5 the state D − A − is only weakly populated and n DA (t) is constant and close to zero. The dynamics is dominated by an electron pair transfer. The combined process is expected if d D eq = d A eq < 0.5. First the population d A rises quickly while n DA stays close to zero. Later a slow increase of n DA to its equilibrium is observed. For the transfer is purely sequential.

IV. THE LOW-TEMPERATURE LIMIT: THE NUMERICAL RENORMALIZATION GROUP
At low temperature, the quantum generalization of the Marcus theory replaces the classical environment by a bath of non-interacting bosonic degrees of freedom. Very early on, the "non-adiabatic" weak coupling limit was investigated. 41 The strong coupling limit of such a model has been adressed using the non-interacting blib approximation (NIBA) 5 , path integral methods 6 and recently also by the numerical renormalization group which we employ in this paper.
Originally the NRG was invented by Wilson for a fermionic bath to solve the Kondo problem. 35,36 The fermionic NRG is a standard and very powerful tool to investigate complex quantum impurity problems. 39 The method was recently extended to treat quantum impurities coupled to a bosonic bath 37,38 , to a combination of fermionic and bosonic baths 42 , and to the calculation of real-time dynamics out of equilibrium. 33,34,43 The non-perturbative NRG approach has been successfully applied to arbitrary electron-bath coupling strengths. 37,38,42,44 A. Equilibrium NRG The numerical renormalization group achieves the separation of energy scales by logarithmic discretization of the energy continuum into intervals Only one single mode of each interval couples directly to the quantum impurity, indicated by the circles in Fig. 3(a). This discrete representation of the continuum is mapped onto a semi-infinite tight-binding chain using an exact unitary transformation. Hereby, the quantum impurity couples only to the very first chain site as depicted in Fig. 3(b). The tight-binding parameters t n linking consecutive sites of the chain m and m + 1 fall off exponentially as t m ∼ Λ −m . Each bosonic chain site is viewed as representative of an energy shell since its energy w m also decreases as w m ∼ Λ −m establishing an energy hierarchy. Both ensures that mode coupling can only occur between neighboring energy shells which is essential for the application of the renormalization group procedure. To this end, the renormalization group trans- . . .
where H m is the Hamiltonian of a finite chain up to the site m -as depicted in Fig. 3 The RG transformation (21) is used to set up and iteratively diagonalize the sequence of Hamiltonians H n . In the first step, only the electronic donor-acceptor system coupling to the single bosonic site m = 0 is considered. It turns out to be sufficient 37,38,39 to include only the N b lowest lying bosonic states, where N b takes typical values of 8 − 12. The reason for that is quite subtle: the coupling between different sites decays exponentially and is restricted to nearest-neighbor coupling by construction, both essential for the RG procedure. In each successive step (i) a finite number of N b bosonic states of the next site m+1 are added, (ii) the Hamiltonian matrices are diagonalized and (iii) only the lowest N s states are retained in each iteration. The discarding of high-energy states is justified by the Boltzmannian form of the equilibrium density operator when simultaneously the temperature is lowered in each iteration step to the order T m ∝ Λ −m w c .
To illustrate the procedure, the lowest-lying energies of the Hamiltonian H m to H m+3 are schematically depicted in panel (c) of Fig. 3. We typically use N b ≥ 8 and keep about N s = 100 states after each iteration using a discretization parameter Λ = 2.
Denoting the set of low-lying eigenstates by |r N and the corresponding eigenvalues E r (N ) ∝ O(1) at iteration N , the equilibrium density matrix ρ 0 is given 39 by where Z N = r e −βE N r andβ are of the order O(1), such that T N = w c Λ −N /β. The thermodynamic expectation value of each local observableÔ is accessible at each temperature T N by the trace Ô eq = Tr ρ 0Ô = 1 Z N r e −βE N r N r|Ô|r N .(23) The procedure described above turns out to be very accurate because the couplings t m between the bosonic sites along the chain are falling off exponentially, so that the rest of the semi-infinite chain contributes only perturbatively 35,39 at each iteration m while contributions from the discarded high-energy states are exponentially suppressed by the Boltzmann factor.

B. Time dependent NRG
While the equilibrium properties are fully determined by the energy spectrum of the Hamiltonian, the nonequilibrium dynamics requires two conditions: the initial condition encoded in the many-body density operator ρ 0 and the Hamiltonian H f which governs its time-evolution. For a time-independent Hamiltonian, the density operator evolves according toρ(t > 0) = e −iH f t ρ 0 e iH f t . All time-dependent expectation values Ô (t) are given by We obtain the density operator ρ 0 from an independent NRG run using a suitable initial Hamiltonian H i . By choosing appropriate parameters in H i , we prepare the system such that (for the calculations presented in this paper) the two electrons are located on the donor site and the acceptor site is empty.
In general, the initial density operator ρ 0 contains states which are most likely superpositions of excited states of H f . For the calculation of the real-time dynamics of electron-transfer reactions it is therefore not sufficient to take into account only the retained states of the Hamiltonian H f obtained from an NRG procedure. The recently developed time-dependent NRG (TD-NRG) 33,34 circumvents this problem by including contributions from all states. It turns out that the set of all discarded states eliminated during the NRG procedure form a complete basis set 33,34 of the Wilson chain which is also an approximate eigenbasis of the Hamiltonian. Using this complete basis, it was shown 33,34 that eq. (24) at iteration m, where ρ 0 is given by (22) using H i . The restricted sum trun r,s in eq. (25) implies that at least one of the states r and s is discarded at iteration m. Excitations involving only kept states contribute at a later iteration and must be excluded from the sum.
As a consequence, all energy shells m contribute to the time evolution: the short time dynamics is governed by the high energy states while the long time behavior is determined by the low lying excitations. Dephasing and dissipation is encoded in the phase factors e i(E m r −E m s )t as well as the reduced density matrix ρ red s,r (m). Discretizing the bath continuum will lead to finite-size oscillations of the real-time dynamics around the continuum solution and deviations of expectation values from the true equilibrium at long time scales. In order to separate the unphysical finite-size oscillations from the true continuum behavior, we average over different bath discretization schemes using Oliveira's z-averaging (for details see Refs. 34,45). We average over 8 different bath discretizations in our calculation.

V. EQUILIBRIUM PROPERTIES
In order to gain a better understanding of the nonequilibrium dynamics presented in Sec. VI, we briefly summarize the equilibrium properties of the model given by Eq. (1). It has been analyzed in Ref. [28], where self-trapping (localization) in the single and two-electron subspace was found.
We start with the phase diagram of the two-site model, as shown in Fig. 4. Only for ε = 0 a quantum phase transition of Kosterlitz-Thouless type separates a localized phase for α > α c from a delocalized phase for α < α c . We plot the phase boundaries between localized and delocalized phases in the α-Ũ -plane, both for single-and two-electron subspaces (grey and black line in Fig. 4, respectively).
For the single-electron subspace, the Coulomb repulsion is irrelevant, and the phase boundary does not depend onŨ. The value of the critical coupling strength, α c , is identical to those of the corresponding spin-boson model. The critical value 5,38 of α c depends on the tunneling rate ∆ and reaches α c = 1 for ∆ → 0.
The phase boundary for the two-electron subspace does depend onŨ, which has drastic consequences for the electron transfer process. Imagine that, by a suitable choice of parameters, the system is placed between the two phase boundaries above the single-electron (grey line) and below the two-electron phase boundary (black line) in the area indicated by I in Fig. 4. Then the system would be in the localized phase in the single-electron subspace. However, one additional second electron immediately places the system in the delocalized phase, and one or even both electrons can be transfered. Similarly, a second electron added to the system in the parameter regime of area II shows the opposite behavior: both electrons get localized although a single electron could be transfered.
Note the different values of the α c 's even forŨ = 0 in the single and the two-electron subspace: the coupling of the donor/acceptor system to the bath induces an effective attractive Coulomb interactionŨ eff = −2αω c between the electrons. On the localized side of the transition, the electron tunneling ∆ is renormalized to zero, so that an electron transfer is clearly absent in this regime. This statement holds only for Ohmic dissipation, on which we focus here; deep in the sub-ohmic regime, coherent oscillations have been recently observed even in the localized phase, see Ref. 44. Figure 5 shows results for the double occupation probability as a function of the electron-bath coupling α for differentŨ calculated with the equilibrium NRG. For the symmetric model considered here, the equilibrium probabilities for the double occupation on donor and acceptor sites are equal: d A eq = d D eq ≡ d eq using the observables defined in Eq. (3). The probability of having two electrons at different sites is given by n DA eq = 1 − 2 d eq . The average double occupancy d eq decreases with increasing effective Coulomb repulsionŨ and increases with increasing α. This can understood in terms of the effective Coulomb interactionŨ eff =Ũ −2αω c , renormalized due to the coupling to the bosonic bath. The delocalization/localization phase transition occurs when d eq → 0.5, as can be seen by comparing Fig. 4 and The double occupancy d eq is calculated analytically for α = 0 and arbitrary ∆ andŨ . For T → 0, the d eq approaches . (27) while in the opposite limit, T → ∞, we obtain d eq → 0.25. The low-temperature limit (27) is included as endpoints of the curves in Fig. 5.
Let us now turn to the temperature dependence of d eq . Figure 6 shows results for temperatures between T = 0.004ω c and T = 0.2ω c for several choices of model parameters. Our calculations imply an independent check of the correct t → ∞ behavior in the next section. Additionally, we can make connection to the high-temperature results of section III. For temperatures T ≪Ũ eff the double occupancy d eq is constant as expected from quantum statistics but deviates drastically from the predictions of the Marcus theory given by Eq. (20). The double occupancy d D eq calculated with the NRG approaches the value 0.5/(1 + eŨ eff /T ) for U eff ≈ T . This result indicates that forŨ eff > T Marcus theory is not applicable while low temperature methods like the NRG are valid.

VI. NON-EQUILIBRIUM DYNAMICS
We employ the time-dependent NRG to evaluate the low-temperature time evolution of the local occupancies using Eq. (25) and investigate the influence of different Coulomb interactionsŨ , single-electron hopping matrix elements ∆, couplings between the electronic system to the bosonic bath α and temperatures T between T = 3 · 10 −8 ω c and T = 0.125ω c . The donor/acceptor sub-system is initially prepared in a state with the two electrons placed on the donor site and evolves according to Hamiltonian (1). We calculate the time-dependent expectation values d D (t) = d D (t), d A (t) = d A (t) and n DA (t) = n DA (t) using Eq. (25). These expectation values are related at any time by the completeness relation d D (t) + d A (t) + n DA (t) = 1. The time evolution of n DA (t) serves as criterion to distinguish between direct two-electron transfer and two consecutive one-electron steps. If n DA (t) remains close to zero or stays constant throughout the electron transfer process, the two states D − A − are only virtually occupied, and concerted twoelectron transfer is observed. A significant increase of n DA (t) as function of time is taken as an indication of step-by-step single-electron transfer.
In the absence of the electron-boson coupling (α = 0), the dynamics is fully determined by the dynamics of the four eigenstates of H el . In the limit |Ũ | ≫ ∆ we obtain while in the limit ofŨ = 0 A finite value of the coupling, α = 0, gives rise to damping of those coherent oscillations. Furthermore, the Coulomb interaction is renormalized toŨ eff =Ũ − 2αω c . ForŨ eff < 0, the states D 2− A and DA 2− are energetically favored. The two intermediate states D − A − are only virtually occupied for |Ũ eff | ≫ ∆, T , similar to the superexchange process. 1 This regime can be described by a spin-boson model with an effective interstate coupling ∆ eff ≈ 4∆ 2 /Ũ eff . The spin-boson model has three dynamical regimes. 5 .
For α smaller than some characteristic value it exhibits damped coherent oscillations between the two states. If α is larger than this value the oscillations disappear and the kinetics is dominated by a relaxation process. Here, rates can be defined and the population probabilities can be fitted with the kinetic equations (19). For a further increase of α, the electronic system shows the onset of localization (for T → 0) and does not evolve towards the other (acceptor) site.
In Fig. 7 we plot the low-temperature population probabilities d D (t), d A (t) and n DA (t) forŨ = −ω c , ∆ = 0.1ω c , T = 3 · 10 −8 ω c and different couplings α. For α = 0 (upper panel) the oscillations have two frequencies (see eq. (28)). The electron pair oscillates from donor to acceptor with the small frequency 4∆ 2 /Ũ whereas the fast oscillations with frequencyŨ characterize the virtual occupation of the high lying states (D − A − ). An increase of α leads to damping of the oscillations (middle panel) and relaxation (lower panel). At about α = 0.3 the electron pair gets self-trapped and the system shows a phase transition to the localized phase at T=0 (see Fig. 4). The configuration D − A − is seen not to be involved in the dynamics as n DA is very small and without ascending slope. Since ∆ ≪Ũ eff the state D − A − cannot be populated as longŨ eff ≫ T .
A more complicated behavior is expected within the four accessible electronic states whenŨ eff ≫ ∆ > 0. In this case the delocalized states D − A − have the lowest energy, and sequential transfer is required to reach the equilibrium state. Pair transfer occurs on a smaller time scale. Thus, a combined pair and sequential transfer on two different time scales governs the dynamics for these parameters.
The four panels in Fig. 8 depict the time evolution of the occupation probabilities d D (t), d A (t) and n DA (t) forŨ = ω c , ∆ = 0.1ω c and four different values of α:  (panel(b)) to purely sequential transfer (panel (c) and panel (d)) with a complex dynamics is due to a combined effect of dissipation and decrease of the effective energy difference between the donor/acceptor states. An even larger α leads to a neg-ativeŨ eff and a very slow transfer until the onset of localization at α c which is not shown here.
To separate the influence of dissipation from the renormalization ofŨ due to the coupling to the bosonic bath we plot d D (t) and n DA (t) for a constant effectivẽ U eff = ω c and different coupling α in Fig. 9. The dynamics changes from pair transfer with a slow increase of the single occupancy at α = 0.04 (due to the low-lying states D − A − ) to incoherent relaxation and sequential transfer for α = 0.36. As long as E α1 = 2αω c ≪Ũ eff , pair transfer is observed on a short time-scale. For E α1 ≥Ũ eff only one electron is transferred and the system relaxes rapidly into its equilibrium state D − A − without any short time concerted pair transfer.
In Fig. 10 the evolution of the dynamics is shown for U = 0 and increasing α. The doubly occupied states are the ground states of the donor/acceptor system for finite α sinceŨ eff = −2αω c < 0. With increasing α, the amplitude of coherent oscillations acquire a small damping. In addition, pair transfer is favored and n DA (t) decreases. The simple damped oscillations are replaced by a much more complex dynamics comprising of strongly renormalized oscillation frequency and a strong damping for α = 0.04. At about α = 0.36 -not shown here -the critical coupling α c is reached and the system is localized.
In the lower damping case (Fig. 11), the transfer is reflected by damped electron pair oscillations forŨ = −ω c in Fig. 11 (a). IncreasingŨ = −0.5ω c in Fig. 11 (b) leads to an increase of the population probability of D − A − and to a change of the fast oscillations whit an approximate frequency ofŨ eff . WhenŨ becomes positiveŨ = 0.5ω c (Fig. 11 (c)) the single-electron transfer becomes fast and the main process unlessŨ eff is not too large. In fact at U = ω c the rate from D 2− A to D − A − becomes smaller ( Fig. 11 (d)) and additional electron pair transfer is observed. The graphs Figs. 11(a) and (d) can be understood in terms of Eq. (28) since |∆/Ũ | ≪ 1 and α = 0.04 is small. By the weak coupling to the environment, ∆ is slightly reduced, and the oscillation amplitude decays exponentially. The difference between the two panels (a) and (d) arises from (i)Ũ eff =Ũ − 2αω c instead of the |Ũ | entering Eq. (28) and (ii) from the dissipation which favors the relaxation into the new thermodynamic ground state: while the oscillation frequencies are roughly the same for |Ũ | = ω c , the delocalized states have a lower energy in Fig. 11(d) so that n DA (t) has to increase to its new equilibrium value. The approximations made in Eq. (28) do not hold any longer for the parameters in Figs. 11(b) and 11(c). The electronic dynamics is governed by additional frequencies and becomes more complex. However, the results can still be analyzed and un- derstood within the analytical results of d D (t), d A (t) and n DA (t) for α → 0.
When the coupling α is increased to α = 0.16, a different picture emerges. Very high frequency oscillations with a small amplitude are superimposed on a slowly decaying d D (t) depicted in in the upper panel Fig. 12. Averaging over those oscillations, we can fit the population probabilities to the kinetic equations (19). By this procedure, we extract the phenomenological rates as function ofŨ eff for fixed α = 0.16. As shown in the lower panel of Fig. 12 the concerted transfer rate k pair [D 2− A→DA 2− ] increases with increasingŨ eff (Ũ eff < 0). This was expected from the rate equation (17) in the classical limit.
The transfer is found to be incoherent and sequential in the higher damped case α = 0.36 for not too largeŨ . The population probabilities are shown forŨ = 0.6ω c , 1.7ω c and 2.5ω c in the upper panel of Fig. 13. By fitting the curves with the help of the kinetic equations, eq. (19), we obtain the rate of the single-electron transfer D 2− A to D − A − which is a non-monotonic function ofŨ eff with a maximum atŨ eff = E α1 ≈ 0.72ω c (see lower panel). It is plotted together with the Marcus rate at T = 0.008ω c (For varying temperatures we found that the fitted rate is approximately constant for temperatures T < 0.008ω c in the considered parameter space.). Although the qualitative behavior is captured by the Marcus rate the asymmetric shape of the NRG result is more realistic in the nuclear tunneling regime. AsŨ increases further the sequential transfer becomes negligible in the inverted region. As a matter of fact, an increasing value ofŨ shifts the system away from the phase transition line deeper into the delocalized phase as can be seen in the equilibrium phase diagram Fig. 4. Here, the dynamics is dominated by coherent pair oscillations with a very small frequency, displayed forŨ = 5ω c in the middle panel of Fig.13. Finally, the effect of temperature is studied in Fig. 14 whereŨ = −0.01ω c , ∆ = 0.001ω c . The temperature is varied from 3 · 10 −8 ω c to 0.125ω c . For T = 3 · 10 −8 ω c to T ≤ 0.02ω c the population probability is temperature independent. As long asŨ eff > T , pair transfer is observed (the probability of D − A − stays constant). As T >Ũ eff the states D − A − are seen to contribute and are thermally populated.

VII. SUMMARY AND CONCLUSION
In this paper, we have studied the electron transfer properties of two excess electrons in a redox system modeled as a dissipative two-site Hubbard model -a model which can be viewed as the simplest generalization of the spin-boson model to include many-particle effects. These many-particle effects are due to on-site and intersite Coulomb interactions, U and V respectively, as well as the effective interactions induced by the coupling to a |Ũ eff | ≤ ∆, single-electron transfer T >Ũ eff single-electron transfer U eff > 0: electron-pair transfer (in addition slow single-electron transfer) U eff < 0:  common bosonic bath. These interaction parameters can be calculated by ab initio methods for a specific system (see, for example, 22,24 ). In our two-site model only the differenceŨ = U − V enters the dynamics. In the presence of a bosonic bath, the effective energyŨ is renormalized toŨ eff =Ũ − 2αω c . An effective attractive inter-actionŨ eff < 0 favors the localization of two electrons on the same site, a repulsiveŨ eff > 0 favors the distribution of electrons on different sites. The intricate correlated dynamics of two electrons depends on the activation energy. This is because, for the tunneling of an electron between two states, energy fluctuations are necessary for the reorganization of the donor-acceptor system and is influenced by an energy difference between the states. Therefore, the transfer characteristics in the unbiased case depends strongly on the effective on-site Coulomb repulsionŨ eff . Three rates have to be considered: the forward and backward rate between the double occupied states (D 2− A, DA 2− ) and the two intermediate degenerate states (D − A − ) as well as the direct rate between D 2− A and DA 2− . How these rates depend ofŨ eff is summarized in Table I. We have performed calculations for the probabilities P (t) of doubly and singly occupied donor and acceptor states using the time-dependent Numerical Renormalization Group method 33,34 . This information helps us to identify conditions under which the systems performs (a) concerted two-electron transfer, (b) uncorrelated, sequential single-electron transfer or (c) fast concerted two electron followed by a single-electron transfer. With the time-dependent NRG method we can describe the crossover from damped coherent oscillations to incoherent relaxation as well as to localization (at T → 0). The temperatures are chosen to be 0.1ω c > T > 3 · 10 −8 ω c . For larger temperatures, when the bosonic bath can be treated classically, the Marcus rates are applicable.
ForŨ eff ≫ ∆, E α1 , T concerted electron transfer occurs in both methods: in the nuclear tunneling regime within the NRG as well as in the limit of a classical bath within the Marcus theory. As long as T <Ũ eff , however, thermal activation is absent and nuclear tunneling is the main process. Only a full quantum mechanical calculation yields the correct relaxation rates which are governed by quantum-fluctuation, dephasing and energy exchange with the environment.
For small ∆/|Ũ eff | we found an effective pair hopping via virtual population of the low lying or high lying states D − A − . When the equilibrium probability for the states D − A − is finite, a slow single-electron accompanies the faster pair transfer. In contrast to the single-electron transfer with a frequency of the order ∆, the frequency of the pair transfer is of the order 4∆ 2 /|Ũ eff |.
The concerted transfer becomes more uncorrelated and sequential at short times at high temperatures (T > U eff ), increasing coupling to the bosonic bath (E α1 ≥ U eff ) or larger single-electron hopping (∆ ≥Ũ eff ). The sequential transfer rate is non-monotonic with increasing U eff . At first, the transition rate from D 2− A to the delocalized states D − A − increases for smallŨ eff > 0, reaches a maximum forŨ eff = E α1 before it decreases again. The rate for the consecutive process D − A − → DA 2− , however, decreases with increasingŨ . For a negative ef-fective Coulomb matrix elementŨ eff , the transfer rate of the second process D − A − → DA 2− is maximal for U eff = −E α1 . In this parameter regime we expect that the second electron follows very shortly after the first electron was transferred.
The transfer kinetics of more than two excess charges in, for example, biochemical reaction schemes or molecular electronics applications is controlled by the molecule specific Coulomb interaction and its polar environment. Our study reveals the conditions for concerted twoelectron transfer and sequential single-electron transfer. Concerted two-electron transfer is expected in compounds where the difference of the inter-site Coulomb repulsion and effective on-site repulsion are much larger than the single-electron hopping and larger than the temperature and reorganization energy. Furthermore, we have shown that non-monotonic characteristics of sequential single-electron transfer strongly depends on the Coulomb interaction. A further study will include the influence of a finite energy difference ε between the donor and acceptor site. We will also report on the influence of Coulomb repulsion and many-particle effects on the long-range charge transfer using a longer Hubbard chain as bridge between donor and acceptor centers.