Correlated charge transfer along molecular chains

The nonequilibrium dynamics of charges along one-dimensional tight-binding lattices subject to a vibronic background is analyzed in the two-particle sector of a dissipative Hubbard model. The role of symmetries for distinguishable fermions with opposite spins and for indistinguishable ones with polarized spins is revealed and is shown to cause transfer dynamics with non-Boltzmann types of asymptotic distributions. Explicit results are obtained within an improved path integral Monte Carlo approach and are compared against an approximate description in terms of rate equations.


Introduction
The theoretical understanding of charge transfer in molecular structures has been subject to intensive research in the last few decades [1]. The most prominent example is the photo-induced primary charge separation in the reaction center of photosynthetic bacteria and plants, where since the pioneering work of Marcus [2] the important role of residual molecular degrees of freedom coupling to the charge transfer has been elucidated in the classical as well as in the quantum regime [3]- [5]. In short, the environment tends to localize the charge and a transfer is energetically forbidden unless environmental fluctuations, associated with rearrangements of e.g. vibron modes, open a channel towards a neighboring site. The details of this simple picture very much depend on temperature though, such that in the low-temperature domain quantum mechanical nonlocality allows for tunneling processes of collective environmental modes and moreover renders, even for modes much faster than the bare charge dynamics, retardation effects to occur. Accordingly, the charge transfer follows a stepwise, i.e. sequential, diffusive dynamics only if its phase coherence is destroyed sufficiently fast.
In the last few years, the subject has gained renewed interest with the advent of new types of electronic devices, where molecular structures are embedded into mesoscopic circuits [6,7]. The nonlinear conductance properties of these devices are intimately related to the charge transfer processes discussed, e.g. in the context of photosynthesis [8]- [10]. A detailed description, however, is much more challenging as it must treat the fermionic nature of the leads with a continuous density of states on equal footing with the electronic and vibronic molecular structure. Artificial molecular structures have been fabricated in the form of coupled quantum dots [11], which allow one via gate voltages to control the occupancy of the individual dots and the transfer between them [12]. While in native structures typically only single charges have been observed, in the described devices correlations between the charges have a substantial impact and may lead e.g. to Pauli-blocking and Coulomb blockade.
Essential properties of environmental-mediated charge transfer have been captured within spin-boson types of models [13], where charges move along a tight-binding lattice interacting with a bosonic heat bath, the spectral density of which mimics the presence of vibrons and/or solvent degrees of freedom. In the challenging low-temperature domain, a particularly elegant approach to describe the nonequilibrium dynamics is the path integral formulation of the reduced density matrix ( [14] and references cited therein). The latter carries information about the charge dynamics after the heat bath has been traced out so that its effective impact appears as an influence functional describing self-interactions of the charges nonlocal in time. This formulation is formally exact and in particular non-perturbative in the coupling between heat bath and charges. Due to the retardation, however, simple equations of motion for the reduced density do not exist and the commonly used master equations are only applicable for sufficiently weak dissipation and high temperatures. Explicit evaluations thus necessitate numerical techniques. In this context, path integral Monte Carlo (PIMC) algorithms have been developed in the last few years [15]- [17] into a powerful means to investigate not only single charge [16,18], but also correlated charge transfer dynamics [19] along finite molecular chains. In this latter situation, the spin-boson model is generalized to a dissipative Hubbard model.
Recently, we have shown [19,20] that symmetries in the bath-charge coupling strongly influence dynamical properties and even lead asymptotically to non-Boltzmann equilibrium distributions. The details of this behavior are determined by the particle statistics. So far the two-charge/opposite spin sector (distinguishable fermions) has been analyzed; to obtain a general picture, however, the situation of spin-polarized charges needs to be considered as well. Moreover, this case reveals the intimate relation between particle statistics and nonequilibrium quantum dynamics in a much deeper way since the Pauli principle drastically reduces the accessible many-body states. Here, we provide this extended account for correlated charge transfer.
The paper is organized as follows. We start in section 2 by briefly recalling the main results for distinguishable charges and introducing the notation. In section 3, the methodology of the PIMC simulations is explained, followed in section 4 by an analysis of the dynamics of indistinguishable charges. As discussed in section 5, in some cases the PIMC data can be well described by rate models based on sequential jumps along the molecular chain. Explicit results are discussed in section 6.

Distinguishable charges
In this section, we briefly recall results for the case of distinguishable charges, which has been discussed in detail recently [19,20]. Accordingly, we concentrate for fermions on the two-particle/opposite spin sector, but mention that all results are fully transferable to any two particles though only as long as they can be physically distinguished. Based on this analysis, the generalizations to more charges/particles may be tedious in detail, but straightforward in principle.
Such a system of two fermions with opposite spin is modeled by an open Hubbard chain with d sites of spacing q 0 according to where c † i,σ and c i,σ are annihilation and creation operators for electrons, respectively, with spin σ on the site i. E i are the bare energies of the levels, U ik denotes the strengths of the Coulomb interaction between sites i and k, and i are the tunneling matrix elements. The influence of the bosonic bath is given by a Caldeira-Leggett type of model [21], where the interaction with the bath is accomplished via a standard dipole coupling [1] where is the polarization operator of the Hubbard chain. In this way, the above model can be seen as a generalization of the spin-boson model [13,14] of charge transfer to the case of many electrons and many sites.
A suitable basis for the electronic degrees of freedom is now given in terms of the localized many-body states (LMBS) where |s ↓ , s ↑ A denotes an antisymmetrized state with the σ = ↓ (σ = ↑) fermion being localized on site with, e.g., |1, ↓, s ↓ denoting a state with the first particle being in the spin-down state, localized on site i = s ↓ + S + 1. For convenience, the subscript A will be omitted in the following. The time evolution of the LMBS populations then follows from where ρ(t) = tr B {W (t)} denotes the reduced density matrix of the electronic system, obtained from the full density matrix W (t) = exp(−iH t/h)W (0) exp(iH t/h) with H = H S + H B + H I after integrating out the environmental degrees of freedom. As shown in [19,20], the reduced dynamics (7) is strongly determined by the symmetries of the underlying Hamiltonian of the total compound. In fact, the dipole type of systembath coupling can provide symmetries which decompose the electronic Hilbert space into disjunct subspaces. The projector onto these subspaces commutes with H S and H I , i.e. with the polarization P, meaning that it is not affected at all by the dissipative environment. There is even one particular subspace corresponding to a vanishing eigenvalue of the polarization. It is thus completely decoupled from the heat bath and is called a decoherence free subspace (DFS). For example, when populating initially only one of these subspaces, the dissipative system's dynamics remains completely restricted to it. This has profound consequences: (i) the equilibrium populations of the invariant subspaces explicitly depend on the initial preparation; in particular, the on-site populations are no longer Boltzmann distributed and (ii) for appropriate bath parameters, the full quantum dynamics can be expressed into simple rate equations, which intrinsically obey the corresponding symmetries and lead to the correct equilibrium populations.
The analysis of the symmetry properties and their relation to invariant subspaces is most conveniently done by exploiting that the dissipative n-particle dynamics on a one-dimensional lattice can be mapped onto an effective single particle diffusion on an n-dimensional one. For instance, as shown in [22], for two distinguishable fermions on a one-dimensional lattice with d sites, a mapping exists onto one-particle dynamics onto a two-dimensional square lattice with d 2 sites, see figure 1. This mapping turns out to be a very convenient tool to visualize the relation between many-body dynamics and dissipation and may further serve as a starting point to apply perturbative methods e.g. the noninteracting blip/cluster approximation (NIBA/NICA), which is known to be a powerful means to capture the single-particle motion. The rate equations mentioned above are a direct consequence of this type of approximation.
For a more in-depth discussion, we first note that, as discussed in [20], the Hamiltonian (1) obeys spin-permutation symmetry due to the fact that no spin-flip or magnetic interactions are included. Consequently, accessible many-body states split into those with even and those with odd eigenvalue with respect to a spin-flip operation. In the case of two fermions, this amounts to the fact that singlet and triplet states do not communicate with each other even in the presence of a bath. Secondly, for the tight-binding lattice with the localized basis (5), it is suggestive to identify the discrete positions of the charges s ↓ and s ↑ with the discrete eigenvalues of the z-components J ↓ z and J ↑ z , respectively, of two pseudo-angular momentum operators. Accordingly, the latter ones obey J ↓2 = J ↑2 =h 2 S(S + 1) with d = 2S + 1. For single charge transfer on two sites this leads to the famous spin-boson model where the polarization is then proportional to the Pauli matrix σ z . Now, for two charges the polarization P is proportional to J ↓ z + J ↑ z . Hence, for the analysis of symmetry properties it is convenient to introduce an alternative basis of the fermionic Hilbert space following from a total pseudo-angular momentum These collective states are delocalized and allow us to represent the original Hamiltonian in a new basis, such that the symmetry encoded in the system-bath coupling is optimally reduced (see figure 2). Accordingly, since [H S , J z ] = 0, invariant subspaces exist for [H S , J 2 ] = 0, each one exhibiting completely isolated dynamics. These subspaces are then identical to the irreducible representations of the corresponding pseudo-angular momentum characterized by its eigenvalue j. While the latter relation requires indeed a specific form of the site-couplings and energies-e.g. if H S is of the form which also allows for finite U -even in the case where [H S , J 2 ] = 0 is not fulfilled, the spin-flip symmetry is always conserved, and subspaces with even j do not couple to those with odd j. Let us discuss this in detail. Within the pseudo-angular momentum representation the localized states in B loc (5) are simultaneous eigenstates of J ↓2 , J ↓ z , J ↑2 and J ↑ z due to while, as already mentioned, the polarization operator (4) can be expressed as Here, J z denotes the z-component of J . An alternative electronic basis set is then given by the simultaneous eigenstates of these two operators, i.e. where Note that, like the localized states in B loc , the collective states in B J are also eigenstates of the polarization operator P (12), as well as of a spin-permutation operator.
The transformation between the two basis sets is given in terms of the (real-valued) Clebsch-Gordan coefficients [23] and exploiting s ↓ , s ↑ |ψ ( j) With the aid of equation (15) one can now easily express the LMBS populations (7) in terms of the corresponding ones in the pseudo-angular momentum basis, yielding Equation (17) also governs the transformation of the initial populations between the two basis sets. Particularly simple expressions are gained for an initially factorizing system-bath state 7 where, as in many experiments, the system is prepared to reside on some localized state |s ↓ 0 , s ↑ 0 . The corresponding initial populationsP j,m (0) are then obtained from For long times the system relaxes towards an equilibrium state. Depending on temperature, spectral density and the system-bath coupling this stationary state may be a thermal equilibrium where all states are occupied according to a Boltzmann distribution. For single charge transfer, this is particularly true for an ohmic spectral density with sufficiently strong coupling and at higher temperatures [14]. For a dissipative Hubbard model, however, the same only holds for each of the disjunct subspaces separately, rather than for the total Hilbert space. In particular, one obtainsP where the initial population in subspace j isP The same is true for the states in B loc . Since in thermal equilibrium all cross terms in equation with Z j = j m =− j e −hβˆ ( j) m . Apparently, the equilibrium populations strongly depend on the initial preparation, i.e. onP

Quantum Monte Carlo simulations
As already addressed above, the path integral formulation of quantum mechanics has been proved to allow for an exact elimination of the reservoir degrees of freedom and to provide, in combination with Monte Carlo techniques, a numerically exact evaluation [3,15]. For single charge transfer along molecular chains, the method has been substantially improved in the last few years to capture also long chains, impurities and external driving [17,18]. Due to the mapping of the many-body transport onto an effective one-particle diffusion on a higher dimensional lattice, also correlated transfer could be simulated. To keep the presentation to some extent self-contained, we collect the main findings in this section.
The expectation value of a system observable A requires the knowledge of the reduced density operator Here, H denotes the full Hamiltonian (1)- (3), W (0) is an initial density, and the trace is performed over the bath degrees of freedom only. For further details of the path integral representation for ρ(t), the reader is referred to the literature [13,14,16,18]. The reduced density operator is expressed as a double path integral along a Kadanoff-Baym contour within a proper basis of the many-body Hilbert space. The straightforward choice uses forward (s σ ) and backward (s σ ) paths corresponding to the basis set B loc . The impact of the dissipative environment appears as an influence functional introducing arbitrarily long-ranged interactions in time along and between the paths. It is convenient to switch to the sum and difference coordinates η σ = s σ +s σ and ξ σ = s σ −s σ , respectively, so that one arrives for the expectation value at the exact expression where a[ η, ξ ] is the measurement functional corresponding to A in terms of the combined system paths η(t) = (η ↑ (t), η ↓ (t)) and ξ (t) = (ξ ↑ (t), ξ ↓ (t)). Furthermore, K is the bare action factor in the absence of a reservoir, and the influence functional reads with e = (1, 1). The kernel is related to the force-force auto-correlation function of the bath and is completely determined by the spectral density J (ω) of its modes [14]. Further, µ = limh β→0h β L(0). Note that even in the absence of Coulomb interaction the two charges are correlated due to the coupling to the common heat bath. An analytical treatment of the expression (24) is in general not feasible, mainly due to the retardations in the influence functional which grow with decreasing temperature, roughly ashβ. In this situation, PIMC methods have been shown to be very powerful and numerically exact means to explore the nonperturbative range. A striking difference from the single-particle case, however, is the existence of symmetries as described above and the consequent decomposition of the full system's Hilbert state into invariant subspaces. To exploit these symmetries for our numerical studies, we may express the path integral equation (24) in terms of system paths J (t) and M J (t), referring to the collective states |ψ (J ) M (13) rather than in terms of the localized ones s ↓ (t) and s ↑ (t). It turns out that the environmental influence can again be summarized in an influence functional which is obtained from equation (25) Note that for fixed J this is just the influence functional of a single dissipative particle residing on 2J + 1 discrete states. Now, in the specific cases that J (t) is conserved during the propagation ([H, J 2 ] = 0), correlations between the invariant subspaces can only be mediated by the initial preparation Tr{W (0)} and the measurement a[η M , ξ M ] at time t. In terms of paths of J and M, equation (24) therefore eventually becomes where J 1 and J 2 denote the values of the piecewise constant J paths along the forward and backward parts of the Kadanoff-Baym contour, respectively, and a J 1 ,J 2 [η M , ξ M ] is the measurement functional of the operator P J 1 AP J 2 with the projector P j onto the invariant subspace j. For example, if the system is initially prepared in the localized state | − S, −S , one simply obtains (29) Accordingly, in these cases the dissipative two-particle system can be decomposed into d independent dissipative single-particle systems. For PIMC simulations, this presents a convenient simplification since each of the systems representing the A J 1 ,J 2 (t) contributions exhibits a significantly lower dimensionality than the full d 2 dimensional two-particle system. Hence, by evaluating equation (28) rather than equation (24), the dynamical sign problem [15], which stems from interferences between different quantum paths and grows exponentially with the size of the system, is correspondingly reduced.
If, due to [H, J 2 ] = 0, a full decomposition into d disjunct subspaces is not possible, often one still retains a partial decomposition where only some of these subspaces are connected to each other. Then, an expression similar to equation (28) can still be obtained, e.g. by exploiting the spin-permutation symmetry. Therefore, the arguments presented above concerning the computational costs of PIMC simulations still apply, albeit to a lower degree.
Technical details of the approach and particularly of the strategy to soothe the dynamical sign problem have been discussed in [16]- [18]. The basic ingredients are as follows: (i) one exploits the fact that the influence functional depends only linearly on the quasi-classical paths η(t), so that in equation (24) the corresponding summations can be expressed as a series of simple matrix multiplications and therefore can be carried out explicitly; (ii) for the sampling procedure one chooses an MC weight with shortened long-time retardations. Since they are fully taken into account in the final accumulation process evaluating the exact expression (24), the numerical exactness of the MC scheme is not impaired. In this way, one achieves a strong decoupling of quasi-classical ( η) and quantum ( ξ ) coordinates, which allows one to store products of short time propagators independent of the MC sampling. Eventually, a speed-up with respect to the original method [16] by a factor of about 100 is gained.

Indistinguishable charges
Now let us turn to situations where charges are indistinguishable, which in the case of fermions means that carriers are spin polarized. In a Hubbard type of Hamiltonian (1), spin degrees of freedom are thus dropped and an on-site Coulomb interaction may be replaced by a nonlocal interaction, i.e.
The mapping from two-particle diffusion along a one-dimensional chain with d sites onto an equivalent one-particle diffusion now identifies all localized states that differ only due to the permutation of particles. Further, any double occupancy is forbidden by the Pauli principle. Thus, the two-body system is mapped onto a triangular two-dimensional lattice with d(d − 1)/2 sites (cf figure 3). We note in passing that formally this coincides with the lattice obtained for indistinguishable bosons diffusing along a chain with d + 1 sites. Instead of the localized basis set (5), we now have where |s, s denotes an antisymmetrized state with one fermion residing on site s and the other on site s with |s, s = −|s , s (cf (6)). Invariant subspaces can still exist. One way of constructing them starts from the representation (15) for distinguishable particles and applies the symmetry property of indistinguishable ones. Then, only those subspaces survive for which Clebsch-Gordan coefficients are antisymmetric with respect to permutations (s, s ) → (s , s), which by virtue of C s ,s means that j − 2S = j − d + 1 (0 j 2S) must be odd. In particular, a DFS, i.e. a subspace with j = 0, only exists if d − 1 is odd. Accordingly, the eigenstates for a subspace with eigenvalue j follow for a chain with d sites from the localized states ( j − 2S odd) as where the prefactor N jm is the corresponding normalization. The lower limit of the sum ensures that only states |s, s with s > s contribute. Similarly to the case of distinguishable charges, the thus constructed subspaces for given j turn out to be disjunctive only if the Hubbard Hamiltonian exhibits well-defined relations between its on-site energies, tunnel couplings and Coulomb terms. A violation of this structure, however, does not necessarily imply the complete absence of any invariant subspaces. For instance, in the case d = 4 a DFS exists according to (32) only if 1 = 3 and is then given by However, for 1 = 3 a DFS can still be found if a/b = 3 / 1 , which is not included in (32).
The Hamiltonian (30) can now be represented in the collective basis to account more efficiently for symmetries when performing PIMC simulations. Since the bath does not induce transitions between invariant subspaces, they again evolve isolated in time. Thus, what has been said in section 3 also applies here.
Further, the LMBS populations P s,s (t) are related to the collective state populationsP j,m according to (17) upon replacing s ↓ , s ↑ by s, s and performing the sum only over those j for which j − 2S is odd. Likewise, the equilibrium populations are given by (19).

Rate description
For single charge transfer along one-dimensional tight-binding systems, a description of the population dynamics in terms of a retarded master equation is completely equivalent to the path integral representation. For sufficiently strong dissipation and/or high temperatures, phase coherence is lost after each jump to an adjacent site and one speaks about sequential transfer. In this case, the exact dynamics obeys approximately where P collects the single-particle populations and the matrix R contains the transition rates between adjacent sites. These rates obey detailed balance reflecting the existence of a thermodynamic equilibrium approached for long times, where the populations are Boltzmann distributed such that R P(t → ∞) = 0. However, as we have shown in the previous sections, for the many-body time evolution this is no longer true, and formulating equation (34) with P(t) = (P −S,−S (t), P −S,−S+1 (t), . . .) and the corresponding golden rule rates in R must inevitably fail to reproduce the correct dynamics. For many-body transfer, one thus starts, in the case of distinguishable charges, from the pseudo-angular momentum basis (15), while for indistinguishable ones, (32) offers an appropriate set for specific types of Hamiltonians. As discussed above, in this latter case invariant subspaces may also exist for other types, where the corresponding many-body states must then be calculated directly from the Hamiltonian. For each invariant subspace, one employs separate rate descriptions with proper initial and equilibrium populations, respectively, and with transition rates chosen according to known expressions like golden rule formulae [17,18]. The LMBS populations can then be obtained from equation (17) and its corresponding form for indistinguishable fermions. This approach, however, only applies for LMBS populations for which the cross terms in equation (17) Figure 4. Initial dynamics of a two-particle system (distinguishable) for d = 3 with completely degenerate on-site energies and coupling strengths 1 = 2 = and U C = 0. Shown are PIMC data for the localized populations equation (7) as obtained from equation (24) (empty black symbols) and equation (28) (filled gray symbols); gray lines denote the results from the rate approach described in section 5. Shown are P −1,−1 (circles), P −1,0 (triangles up), P −1,1 (triangles left), P 0,−1 (triangles down), P 0,0 (squares), P 0,1, ('+'), P 1,−1 (triangles right), P 1,0 ('×') and P 1,1 (diamonds), with equilibrium values of 1/5, 1/10, 1/30, 1/10, 2/15, 1/10, 1/30, 1/10 and 1/5, respectively.
m } (indistinguishable), respectively. In terms of the path-integral picture, this restricts a rate approach to populations whose path-integral expression (28) collapses to a single A J 1 ,J 1 (t).
For instance, for a distinguishable system obeying [H, J 2 ] = 0 and being prepared according to P −S,−S (0) = 1 or P S,S (0) = 1, this holds for all LMBS populations, while for P 0,0 (0) = 1 only some of them can be reproduced this way (cf figure 4) (for details, see [20]). If the initial preparation does extend over more than one subspace, the dynamics can by no means be captured by a simple rate description along the lines of equation (34). Accordingly, correlations connect disjunctive subspaces and corresponding coherences may survive for longer times, an effect clearly beyond the scope of simple rate models. Thus, Monte Carlo simulations are of particular importance to obtain insight into the nonequilibrium dynamics of correlated charge transfer in the generic cases where symmetries are broken.

Discussion
In this section, we apply the general results of the previous sections to tight-binding systems with d = 3 and d = 4 sites. The simplest case d = 2 is not explicitly addressed since for distinguishable charges it has been extensively studied in [19] and for indistinguishable charges it trivially collapses. We will present only a few typical results for distinguishable fermions and mostly concentrate on spin-polarized charges since the former case has been extensively studied in [20]. In all cases, we worked with a generic spectral bath density of ohmic form with with parameters α = 0.25, ω c = 5 andhβ = 0.1 .

Distinguishable charges
According to the above discussion of isotropic tunnel couplings and all on-site energies including the Coulomb energy set to zero, the full Hilbert space decomposes completely. The rate model adapted to the existence of invariances does indeed capture the dynamics even quantitatively quite accurately. Turning on the Coulomb interaction, i.e. U = 0, a coupling between |ψ (0) 0 and |ψ (2) 0 is introduced (cf figure 2), with ψ (2) 0 |H S |ψ (0) 0 = ( √ 2/3)U , such that only two invariant subspaces remain. However, when fixing the interaction strength to the specific value U C = −1 + 1 − 2 0 , a full decomposition is recovered, see figure 5.

Indistinguishable charges
The number of physically relevant states shrinks significantly when one considers the spinpolarized transfer of two charges. For a chain with d = 3 sites, one recovers the situation of a simple one-dimensional chain with three sites corresponding to the LMBSs |0, −1 , |1, −1 and |1, 0 . The on-site energy of a state |s, s along this effective chain is the sum of the total on-site energies E s + E s and the nonlocal Coulomb interaction U ss . For illustrations, the reader is referred to previous studies for single charge transfer along one-dimensional chains [18]. More interesting is a chain with d = 4 (S = 3/2) sites, which according to our above discussion supports invariant subspaces, in particular a DFS. Let us start with the simplest case E i = U ik = 0, i = and an initial preparation restricted to the subspace j = 2 withP 2,−2 (0) = P −3/2,−1/2 (0) = 1 (see figure 3). The two-body transfer along the chain in subspace j = 2 consists of consecutive jumps of the charges, where only the state |  Accordingly, the populations P s,s (t) for s + s = 0 are given byP 2,s+s (t) and for s + s = 0 one has P 3/2,−3/2 = P 1/2,−1/2 =P 2,0 /2 sinceP 0,0 = 0. From figure 6 one observes that the modified rate model gives for the localized state populations indeed a very accurate description. In particular, the dynamics of P 3/2,−3/2 = P 1/2,−1/2 captures two entangled processes constrained by the Pauli principle, namely, one in which one charge remains on the initial site −3/2, while the other one proceeds to 3/2, and another one in which one charge has already arrived on 1/2 and stays there, while the other one leaves the initial site −3/2 to follow onto site −1/2. The populationP 2,0 shows deviations that can be attributed to the shorter dwell time on this site associated with an incomplete loss of phase coherence due to the fact that this site has a by a factor of √ 2 larger coupling to its adjacent sites. For nonvanishing E s and/or nonvanishing Coulomb interaction U ss , symmetries may be broken and subspaces with j = 2 and j = 0 may be coupled. Thus, populations P s,s with s + s = 0 are determined also by off-diagonal elements of the density matrix according to (17). In general, the modified rate model cannot capture these coherences so that corresponding populations are not well described. In addition, the dwell time on the collective state | (2) 0 is reduced due to the fact that it is linked to three other sites, where the tunnel coupling to the sites within the j = 2 subspace is by a factor of √ 2 larger. The loss of phase coherence may therefore be incomplete. In figure 7, corresponding results for nearest-neighbor Coulomb interaction U i,i+1 = 2.5 verify these expectations. For large Coulomb interaction (U i,i+1 = 40, see figure 8) something notable occurs. Charges initially prepared on adjacent sites are basically trapped and build a type of metastable state stabilized by the phonon background. A similar type of trapping has been also found for distinguishable charges and on-site Coulomb interaction. Namely, as already known from single charge transfer along onedimensional chains, different energy offsets between sites introduce phonon activation barriers between these sites, which is a unique signature of the interaction between the system and phonon bath. This reflects the fact that for a transition between adjacent sites, energy fluctuations are necessary for the reorganization of the bath degrees of freedom after the tunneling of an electron. Here, on the collective state, lattice energy offsets are caused by the Coulomb repulsion that eventually dominates the dynamics such that a fermionic dimer or two-body polaron is formed. We note that, in principle, the methods described here can also be used to investigate the existence and consequences of invariant subspaces for systems consisting of more than two particles. With respect to the PIMC simulations, though, allowing for additional particles will lead to a severly enhanced sign problem, even when exploiting the decomposition of the Hilbert space into disjunct subspaces. However, since the technique used here to battle the sign problem can be regarded as a blocking scheme [18], introducing additional levels of blocking in the spirit of the multilevel blocking approach for real-time calculations [24] can be used to further soothe the sign problem. Finally, the methodology developed here for indistinguishable fermions can directly be used for the dynamics of bosons in a dissipative environment (see the remark above equation (31)), e.g. for the diffusive motion of cold atoms in optical lattices.