Influence of static disorder and polaronic band formation on interfacial electron transfer in organic photovoltaic devices

TL;DR

This study uses quantum dynamical simulations to explore how static disorder and polaronic band formation affect interfacial electron transfer in organic photovoltaics, revealing that dynamic potentials hinder charge separation efficiency.

Contribution

It provides microscopic evidence of how static and dynamic potentials influence charge-injection rates, highlighting the detrimental role of polaronic band formation caused by electron-vibration interactions.

Findings

Dynamic potentials can form polaronic bands.

Dynamic disorder significantly reduces charge-separation efficiency.

Static disorder also impacts interfacial charge transfer.

Abstract

Understanding the interfacial charge-separation mechanism in organic photovoltaics requires, due to its high level of complexity, bridging between chemistry and physics. To elucidate the charge separation mechanism, we present a fully quantum dynamical simulation of a generic one-dimensional Hamiltonian, which physical parameters model prototypical PCBM or $\text{C}_{60}$ acceptor systems. We then provide microscopic evidence of the influence random static and dynamic potentials have on the interfacial charge-injection rate. In particular, we unveil that dynamic potentials, due to strong electron-vibration interactions, can lead to the formation of polaronic bands. Such dynamical potentials, when compared to random static potentials, can provide the main detrimental influence on the efficiency of the process of interfacial charge-separation.

Equations

$$H$$

$$G_{ii}(z{)}^{-1}=z-{ϵ}_i-{\Delta }_i(z)-{\Sigma }_i(z)$$

$$n_{\gamma }(t)=\text{absolutevalue}{\text{innerproduct}{\mathcal{K}}_{\gamma }\psi (t)}^2$$

$$k_{Inj}({ϵ}_0)=-2m^2\text{imaginary}{G}_{11}({ϵ}_0+i\eta )$$

$$k_{Inj}^{\prime }({ϵ}_0)={\mu }^2/\lambda \cdot k_{Inj}({ϵ}_0/\lambda )$$

$$Y\simeq \frac{k_{Inj}}{k_{Inj}+k_{Rec}}.$$

$$K{S}_N(|{\varphi }_0⟩)=span\left\{|{\varphi }_0⟩,H|{\varphi }_0⟩,...,H^{N-1}|{\varphi }_0⟩\right\}.$$

$$|{\mathcal{K}}_0⟩,|{\mathcal{K}}_1⟩,|{\mathcal{K}}_2⟩,\dots ,|{\mathcal{K}}_{N-1}⟩\equiv K{S}_N$$

$$H|{\mathcal{K}}_n⟩=a(n)|{\mathcal{K}}_n⟩+b(n)|{\mathcal{K}}_{n+1}⟩+b(n-1)|{\mathcal{K}}_{n-1}⟩$$

$$H_{\mathcal{K}{\mathcal{S}}_{\mathcal{N}}}=\left[\begin{array}{cccccc}a(0)& b(0)& 0& 0& \cdots & \\ b(0)& a(1)& b(1)& 0& 0& \cdots \\ 0& b(1)& a(2)& b(2)& 0& \cdots \\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots \end{array}\right].$$

$$|\psi (t)⟩$$

$$H_{IDMFT}|0,0⟩$$

$$H_{IDMFT}|{\mathcal{K}}_1⟩$$

$$|{\mathcal{K}}_2⟩=\frac{t|2,0⟩}{b(1)}+\frac{b_{{\Sigma }_0}(0)|1,1⟩}{b(1)}.$$

Full text

I Introduction

Organic photovoltaic devices have currently attracted intense interest due to the potentially cheap production and their ability to environmentally friendly power generation li2012polymer ; gunes2007conjugated ; lu2015recent . The efficiency of these organic photovoltaic devices (OPVs) strongly depends on the charge separation process between materials that transport electrons (usually a fullerene derivative) or holes (usually a polymer). The energy offset caused by this interface drives electrons from the donor to the acceptor, while leaving holes left behind. Once the charges are in separate phases, they need to overcome their mutual Coulomb attraction. This binding energy is in the range of $0.1\text{\,}\mathrm{eV}0.5\text{\,}\mathrm{eV}$ zhu2009charge ; drori2008below ; hallermann2008charge , which is much larger than the thermal energy of about $25\text{\,}\mathrm{meV}$ at room temperature. Surprisingly, the charge separation process and eventually the formation of free charges that can be extracted at the electrodes is still very efficient shirakawa1977synthesis ; morel1978high . In particular, experimental studies have shown, that the charge separation occurs on ultra-fast time scales in the $10\text{\,}\mathrm{fs}100\text{\,}\mathrm{fs}$ range grancini2013hot ; jailaubekov2013hot ; gelinas2013ultrafast . At present, the mechanism responsible for charge separation is not well understood and still actively debated.

An obvious question in organic materials relates to the relative importance of static and dynamic potentials and their effect on the charge separation process. The static (time-independent) potential reflects the electron-hole Coulomb interaction, and the spatial disorder caused by electrostatic interactions resulting from the different environment in which each molecule is placed. This is particularly true for organic devices made up of two different disordered materials ballantyne2010understanding ; hoppe2004organic where the spatial disorder is usually originating from the rather large permanent electric dipole of PCBM. However, we stress here, that contrary to PCBM, molecular dynamic simulations have shown, that \ceC60, which has no permanent dipole, exhibits extremely limited disorder which thus can to a good extent be neglected d2016charges ; tummala2015static . On the contrary, the dynamic potential is related to electron-vibration interactions that result in a time-dependent variation of microscopic transport parameters. Though debated gao2014charge ; de2017vibronic , there are indications that, for charge separation, this effect could contribute falke2014coherent ; song2014vibrational . Thus, the consideration of static and dynamic potentials seems unavoidable for a fully microscopic understanding of the charge separation mechanism in OPV devices. However, such microscopic descriptions are computationally challenging. To tackle these challenges, various numerical methods such as exact diagonalization wellein1997polaron ; wellein1998self , diagrammatic Monte Carlo (DMC) de1983numerical ; kornilovitch1998continuous , time-dependent density functional theory (TTDFT) rozzi2013quantum , and other approaches jeckelmann1998density ; barivsic2002variational have been proposed. However, these methods are either quite expensive from a computational point of view or their application to non-translationally invariant systems remains unclear.

To overcome these difficulties the inhomogeneous version of the dynamical mean-field theory approximation (I-DMFT) potthoff1999surface ; potthoff1999metallic ; potthoff1999dynamical ; delange2016large ; backes2016hubbard ; jacob2010dynamical ; turkowski2012dynamical , a powerful non-perturbative technique for strongly interacting systems, has been introduced. By applying the I-DMFT approximation to a generic one-dimensional model Hamiltonian, which parameters model prototypical PCBM and \ceC60 based acceptor systems, we provide a fully quantum dynamical simulation of the charge separation process taking static (disorder + electron-hole Coulomb interaction) and dynamic potentials (electron-vibration interaction) into consideration. This provides the possibility to compute the charge injection rate at the donor-acceptor interface. Our work provides a step forward to a long-standing challenge in OPV, thereby bridging between chemistry and physics. In particular we unveil here, that dynamic potentials (related to polaron formation), when compared to random static potentials, present the main detrimental lose mechanism in OPVs devices. Yet it can in some instances lead to an enhancement of the charge transfer process.

This paper is organized as follows. Section II is divided into two parts: (A) we introduce the generic one-dimensional Holstein based Hamiltonian to model the charge carrier dynamics across organic model interfaces (B) we briefly introduce the used I-DMFT approximation and comment on its general numerical aspects. In section III we apply the I-DMFT approximation to study the charge carrier dynamics across donor-acceptor model interfaces. Finally, we provide a brief conclusion and outlook in section IV.

II Methodology

Inspired on the Holstein model used in bera2015impact , we use the following generic one-dimensional Hamiltonian to describe the microscopic charge transfer process of an electron at the molecular donor-acceptor interface:

[TABLE]

where $c_{i}^{+}$ ($c_{i}$) is the creation (destruction) operator of electrons, $a_{i}^{+}$ ($a_{i}$) is the creation (destruction) operator of phonons, $\Omega$ is the relevant phonon frequency ($\hbar=1$ throughout this work), $g$ is the electron-phonon coupling strength, $m$ is the tunneling amplitude between the LUMO of the donor (site i = 0) and the first site of the LUMO of the acceptor, $V$ determines the Coulomb potential, $J$ is the electron hopping parameter, $\epsilon_{0}$ is the energy of the incoming electron and $\epsilon_{i}$ is the energy level of a molecule on the acceptor site which is taken to be a random variable (with a mean $\epsilon$ and a standard derivation $\sigma_{\epsilon}$) drawn from a rectangular distribution. In this model, the static potential is given explicitly by the distribution of the onsite energies $\epsilon_{i}$ and the Coulomb potential $V$ between the electron and hole, while the dynamic potential is caused by the electron-phonon coupling term $\sum_{i}g\,c_{i}^{+}c_{i}(a_{i}^{+}+a_{i})$.

At this point, we want to note, that we make the following assumptions. First, we do not model the hole dynamics since the effects of hole diffusion usually lead to a reduction of the Coulomb interaction only resulting in an increase of the electron-hole separation yield athanasopoulos2017efficient ; aram2016charge and since hole transfer typically occurs on time scales in the $1\text{\,}\mathrm{ps}2\text{\,}\mathrm{ps}$ range which is several orders of magnitude slower than electron transfer long2014asymmetry . Second, we use a one-dimensional model to describe a complex three dimensional bulk molecular heterostructure, since the main features of the Holstein polaron do not depend strongly on the dimensionality of the system ku2002dimensionality and since we expect, that a realistic, three-dimensional nature of the system will lead to mainly quantitative changes in the quantum efficiency. Third, we assume a single high-frequency intramolecular mode of vibration (order of $1600\text{\,}{\mathrm{cm}}^{-1}$, i.e. a period of roughly $20\text{\,}\mathrm{fs}$) only, although the involvement of multiple phonon modes provides additional transfer channels. However, in this work, we are interested in electron dynamics that concern the fast electron unbinding from the Coulomb well which occurs on time scales of $10\text{\,}\mathrm{fs}100\text{\,}\mathrm{fs}$. Therefore, the impact of weakly coupled low-frequency modes of vibration (order of $\lesssim$80\text{,}{\mathrm{cm}}^{-1}$$, i.e. a period of roughly $415\text{\,}\mathrm{fs}$) can be neglected, since final equilibration occurs at longer relaxation times xie2015full . Fourth, all computations are presented at zero temperature taking no entropic effects into consideration, since electronic and vibrational energy scales are much larger than $k_{\mathrm{B}}T$, where $k_{\mathrm{B}}$ is Boltzmann’s constant and $T$ is the temperature pelzer2016charge , and since we expect that a change in entropy plays a diminished role in the charge separation process in one-dimensional systems gregg2011entropy . The driving force for the charge separation mechanism is thus entirely of quantum mechanical nature steaming from the coupling of an initially discrete state (electron at the interface) to a final state in the continuum (electron on the acceptor side surrounded by a cloud of phonon excitations). Fifth, we note that in absence of electron-phonon interaction scattering on defects can lead to Anderson localization anderson1958absence and it will hinder electron transfer across the interface. Since it is well known, that in the case of one-dimensional infinite disordered systems, any amount of disorder produces Anderson localization, we have embedded the system into an effective medium which has been computed using the coherent potential approximation (CPA). The absence of Anderson localization within CPA Haydock_1974 will then mimic a system without localization at moderate disorder (except very close to the band-edges) throughout this work.

The approach that we propose is the use of the single-polaron I-DMFT approximation. The aim of I-DMFT is to fully address the relevant spatial variations of the physical properties, but still affording a good description of the physical processes of interest. A crucial aspect of I-DMFT is that it provides an interpolation between the non-interacting case, in which it gives the exact solution of the problem and the strong coupling limit in which is also becomes exact. In I-DMFT the local Green’s function is for one given realization of disorder and is defined as:

[TABLE]

where $\Delta_{i}(z)$, $\Sigma_{i}(z)$ are hybridization function, respectively self-energy at site i and $z=E+i\eta$ with $\eta$ being an infinitesimal small number. In standard I-DMFT procedures a self-consistent solutions is obtained iteratively at each z where one needs to compute repeatedly the diagonal of the inverse of a complex matrix, which dimension equals the number of lattice sites. Using conventional linear-algebra algorithms, this problem grows cubic with the size of the system freericks2016generalized .

We use an alternative approach to solve the I-DMFT self-consistency equations which is based on Haydock’s recursion scheme applied to suitably defined Hamiltonians and, in particular, which does not require the inverse of a complex matrix. Instead the Hamiltonian (1) is solved on the full lattice under the approximation that the electron-phonon self-energy is local depending on the frequency of the local phonon only. Self-consistency equations are than expressed in Hilbert-space such that the recursion technique by Haydock haydock1980solid can be used which makes this method immediately generalizable to any lattice geometry and/or disorder distribution while easily handling inhomogeneous systems of up to $10^{3}$ lattice sites in less than $24\text{\,}\mathrm{h}$ (sequential computation on an Intel Xeon E5-2670 processor).

An detailed explanation of the used I-DMFT approach can be found in the original work by richler2018inhomogeneous .

Finally, let us comment on general numerical aspects of the proposed I-DMFT formalism. The number of possible phonon configurations is infinite, but can be restricted to a finite, sufficiently large number in actual calculations by choosing $M>>g^{2}/\Omega^{2}=\alpha^{2}$, where $M$ is the maximum number of phonon excitations per site. Further, we simulate only a finite part of the total size of the system ($N$). We then embedded the leftmost site into an effective medium which has been computed using the coherent potential approximation (CPA). We thus define a number $N_{i}<N$ of non-equivalent lattice sites along the surface. Choosing $N=2000$, $M=60$ and $N_{i}=300$ ensures that all results are independent of any system size characteristics, while keeping modest computational complexity. A pictorial representation of the proposed model is depicted in Fig. 1.

III Results

In the following we express all energies in units of $J$ (energy unit in $J\approx$0.2\text{,}\mathrm{eV} and time unit in $1/J\approx$3.3\text{\,}\mathrm{fs}) and take $m=0.5$, $V=1.5$, $\Omega=1.0$, $\epsilon=0.0$ and $\sigma_{\epsilon}=0.0,0.1$ while keeping $g$ and $\epsilon_{0}$ as independent parameters. All are chosen such, that they are in a realistic experimental range to model the prototypical PCMB ($\sigma_{\epsilon}=0.1$) and \ceC60 ($\sigma_{\epsilon}=0.0$) acceptor systems zheng2017charge ; d2016electrostatic ; antropov1993phonons ; faber2011electron ; castet2014charge . We have checked that all results are qualitatively insensitive to different disorder configurations, a result that has been found in all tested cases throughout this work. Finally we note, that throughout this paper our initial state at time $t=0$ will consist of an electron at site $i=0$ having energy $\epsilon_{0}$ and no phonon modes excited, i.e. $\ket{\psi(t=0)}=c_{0}^{+}\ket{0}$ with $\ket{0}$ being the vacuum state for phonons and electrons.

First, we present in Fig. 2 the probability density of the wave-function $\ket{\mathcal{K}_{\gamma}}$ on the acceptor

[TABLE]

where $\ket{\mathcal{K}_{\gamma}}$ is the orthonormal basis vector of the $N$-dimensional Krylov-space ($KS_{N}$) gutknecht2007brief . The Krylov-states $\ket{\mathcal{K}_{\gamma}}$, which are computed by a Lanczos-based recursion method haydock1980solid , represent a basis of excitations of the many-body system (electron and phonon modes) that progressively spread away from the interface into the acceptor. The time evolution of the wave function in Krylov space is then defined as $\ket{\psi(t)}=e^{-i\mathcal{H}_{KS_{N}}t}\ket{\mathcal{K}_{0}}$ with $\mathcal{H}_{KS_{N}}$ being the reduced Hamiltonian of the Krylov subspace $KS_{N}$. The time-evolution has then been determined by an exact diagonalization technique by choosing the first Krylov-space vector equal to the initial state, i.e. $\ket{\mathcal{K}_{0}}=c_{0}^{+}\ket{0}$ (see Appendix A for a brief introduction to the Krylov subspace method and computation of wave-function $\ket{\mathcal{K}_{\gamma}}$). This provides an efficient way to extract the essential character of the Hamiltonian (1) while using a limited number of basis sets. We stress here, that the size of the system (2000 states in the Krylov-space) is sufficient to prevent the wave packet to bounce back at the boundary.

We find, in the limit of sole dynamic potentials ($\alpha^{2}=0.5$, $\sigma_{\epsilon}=0.0$) shown in Fig. 2(a), that the probability weight near the interface is progressively decaying with time and the particle is fully delocalized. Upon increasing the electron-phonon interaction ($\alpha^{2}=1.5$, $\sigma_{\epsilon}=0.0$), the local density of states (LDOS) fragments into polaronic sub-bands (self-trapping) and where the strongly renormalized width of the polaronic sub-bands arises as the new energy scale ciuchi1997dynamical . As can be seen in Fig. 2(b) self-trapping of the electron hinders, in this case, the interfacial electron transfer drastically as parts of the wave-function remain localized at the interface resulting in a poor but finite interfacial charge-transfer efficiency. In panel (c), (d) of Fig. 2 we now present the combined effect of static spatial disorder and weak ($\alpha^{2}=0.5$, $\sigma_{\epsilon}=0.1$), respectively strong ($\alpha^{2}=1.5$, $\sigma_{\epsilon}=0.1$) dynamic potentials. As can be seen, the in panel (c), (d) presented physical picture does not change quantitatively when compared to panel (a), respectively (b), an outcome we have found throughout this work.

Before proceeding we stress here, that the LDOS is a key factor since its shape and, in particular, the energy distribution of electronic states are found to determine the value of practically achievable injection energies. This comes in handy when using the proposed formalism richler2018inhomogeneous . In particular, we have found, that taking the incoming electron energy outside the band of delocalized states completely suppresses, as expected, interfacial charge transfer (electron becomes localized at interface) while an incoming electron energy that is taken in the band of delocalized states results in interfacial charge transfer.

By making use of the by I-DMFT computed Green’s function $G_{11}(z)$ and by applying Fermi’s Golden rule we now quantify, as a function of the incoming electron energy $\epsilon_{0}$, the quantum mechanical donor-acceptor interface transfer rate

[TABLE]

where $\imaginary G_{11}(z)$ is the imaginary part of the Green’s function $G_{11}(z)$. This present a simple measure of the efficiency for charge transfer at the donor-acceptor interface, but note, the transfer rate is related, but not equal to the internal quantum efficiency which is commonly extracted from experiments. The interface transfer rate is graphically depicted in Fig. 3 for incoming electron energies of $-2J\leq\epsilon_{0}\leq 2J$, where this restriction is sufficient to assure the applicability of Eq. (4) (weak coupling to a quasi-continuum of states). Although, we do not model the competing electron-hole and/or geminate recombination processes, we have indicated in the same figure the range of experimentally observed recombination rate $k_{Rec}$ of order $0.01\text{\,}{\mathrm{fs}}^{-1}0.05\text{\,}{\mathrm{fs}}^{-1}$ grancini2013hot ; jailaubekov2013hot ; gelinas2013ultrafast . Therefore, efficient charge transfer from the donor to the acceptor occurs when the interface transfer rate $k_{Inj}$ is large compared to the competing process of charge recombination $k_{Rec}$. We start considering the case of sole, but weak dynamic potentials only ($\alpha^{2}=0.5$, $\sigma_{\epsilon}=0.0$). We find, that the interfacial transfer rate is essentially, with the exception of values of $\epsilon_{0}$ close to the low energy edge, higher than the upper experimental recombination rate for all negative values of $\epsilon_{0}$. On the contrary, this changes drastically for all positive values of $\epsilon_{0}$, where interfacial charge transfer is large only, when compared to extremely low experimental ranges of recombination rates. This stems from strong relaxation processes (phonon emission) caused by molecular vibrations which are facilitated at higher energies. Upon increasing the electron-phonon interaction ($\alpha^{2}=1.5$) multiple gaps throughout the spectrum arise and thus multiple values of $k_{Inj}$ occur that are below the experimental value of $k_{Rec}$. This is readily explained since upon increasing the electron-phonon interaction the LDOS fragments into sub-bands that are separated by multiples of $\Omega$. This effect will drastically hinder the range of suitable incoming electron energies, since the incoming electron energy $\epsilon_{0}$ must be identical to that of the allowed unoccupied state of the acceptor material. Surprisingly, this picture does not change qualitatively when considering the combined effect of static disorder and dynamic potentials. In particular, we have found that weak but static disorder does not change drastically the interfacial transfer rate for incoming electron energies of $-2J\leq\epsilon_{0}\leq 2J$. Finally, we note that these results are in agreement with the analysis of the time-evolution of the electron density shown above.

We note here, that we confine our analysis about the interfacial transfer rate by choosing realistic experimental values for all parameters. However, introducing $m\to\mu m$, $J\to\lambda J$ enables the scale transformation of the interfacial transfer rate $k_{Inj}\to k_{Inj}^{\prime}$ into

[TABLE]

where $\mu$, $\lambda$ are the scaling constants of $m$, $J$ respectively. Through this scale transformation the results presented in Fig. 3 can then be adapted to describe physically equivalent systems having a different scale $(k_{Inj}^{\prime})$ and set of parameters ($\Omega^{\prime}=\lambda\Omega,V^{\prime}=\lambda V,g^{\prime}=\lambda g,\sigma_{\epsilon}^{\prime}=\lambda\sigma_{\epsilon}$).

Finally, we remark the following: First, we have not found interfacial electron-transfer in the limit of extremely strong electron-phonon coupling ($\alpha^{2}\gtrsim 8.0$). This is readily explained since in this limit the coupling between the donor and the first site of the acceptor becomes much larger than the renormalized bandwidth of the polaronic sub-bands. Thus, one is, to a good extent, left with two eigenstates separated by a larger energy offset giving rise to charge localization at the interface.

Second, we note that the results of our time-dependent study present some similarity with the Dirac-Frenkel time-dependent simulations in one-dimension bera2015impact where their study has focused on the interplay between electron-vibration interaction and the interfacial Coulomb interaction between the hole and electron. Yet, bera2015impact has not taken spatial disorder into consideration and has focused on the quantum yield in absence of any electron-hole recombination. The quantum yield, in absence of recombination, then simply measures the charge injected at infinite distance from the interface, i.e. their quantum yield is simply $1-P$ with $P$ being the weight of the wave-function localized near the interface. Instead, our study quantifies the charge injection-rate which is of more fundamental relevance since this quantity can be compared to the experimentally observed recombination rate. The quantum yield, in the presence of recombination, can then be obtained by

[TABLE]

A high quantum yield is then given when $k_{Inj}>>k_{Rec}$ holds. Moreover, the applicability of the in Ref. bera2015impact presented simulation is limited since their approximated treatment of the many-body nature of the polaronic state (momentum average approximation) is valid in the small bandwidth limit only. Contrary to Ref.bera2015impact , the proposed I-DMFT approach is accurate over the entire parameter space richler2018inhomogeneous .

IV Conclusion

To summarize, we have applied the I-DMFT approximation to a generic one-dimensional model Hamiltonian, which parameters model the charge carrier dynamics in prototypical PCBM and \ceC60 acceptor systems. Our results show that polaronic bands, when compared to spatial disorder, can provide here the main detrimental influence on the efficiency of charge transfer of electrons across organic interfaces. From this perspective, organic molecules with moderate reorganization energies should be used preferentially in next-generation materials since increasing the electron-phonon interaction hinders the range of suitable incoming electron energies due to the fragmentation of the local density of states into narrow polaronic sub-bands.

Finally, we emphasize here that the easy numerical implementation of the I-DMFT approximation richler2018inhomogeneous allows one to study a variety of recently proposed perhaps more realistic donor-acceptor model systems. In particular, I-DMFT enables to investigate the impact electric fields induced by energy level pinning liu2008control , structural heterogeneity as a function of distance to the interface poelking2015design ; mcmahon2011holes and gradients in the energy landscape poelking2015impact ; wilke2011electric . These problems were previously difficult to access but may assist the charge separation process drastically. Open questions we reserve for future work.

V Acknowledgments

The authors thank S. Fratini, S. Ciuchi and G. D’Avino for the stimulating discussions. K.-D. Richler acknowledges the LANEF framework (ANR-10-LABX-51-01) for its support with mutualized infrastructure.

Appendix A Krylov subspace method

We recall here a general definition of the in this work used Krylov subspace method. The Krylov subspace of dimension $N$ ($KS_{N}$) gutknecht2007brief is the linear subspace spanned by:

[TABLE]

Here $\ket{\phi_{0}}$ represent a suitably chosen vector of the original Hilbert space and where we will denote in the following the orthonormalized basis vectors in $KS_{N}$ by

[TABLE]

with $\ket{\mathcal{K}_{0}}=\ket{\phi_{0}}$ being defined as the initially chosen wave-function. An orthogonal basis of $KS_{N}$ can then be constructed with the method of Haydock Lanczos1950zz , an iterative procedure that is capable of constructing the Krylov space via a three-term recurrence relation haydock1980solid :

[TABLE]

with the initial conditions $\ket{\mathcal{K}_{0}}=0$, $b(-1)=0$ and where $\ket{\mathcal{K}_{n}}$ obeys the orthogonality relation $\bra{\mathcal{K}_{n}}\ket{\mathcal{K}_{m}}=\delta_{n,m}$. The reduced Hamiltonian matrix $H_{\mathcal{KS_{N}}}$ in $KS_{N}$ then reads

[TABLE]

The time evolution of the wave function $\ket{\mathcal{K}_{0}}$ in the Krylov subspace $KS_{N}$ is then defined by

[TABLE]

Here $E_{\mathcal{KS_{N}}}(l)$, $\ket{Z_{\mathcal{KS_{N}}}(l)}$ represent the l-th eigenvalue, respectively eigenvector of the Hamiltonian $H_{\mathcal{KS_{N}}}$ which have been determined by exact diagonalization (due to greatly reduced size of the Hamiltonian matrix $H_{\mathcal{KS_{N}}}$).

To further illustrate the extension of the Krylov subspace vector $\ket{\mathcal{K}_{N}}$, we have added in Fig.4 the tight-binding representation of the Hamiltonian $H_{IDMFT}$, which has been computed with the in Ref.richler2018inhomogeneous presented I-DMFT method, and which has been used to determine all $\ket{\mathcal{K}_{N}}$. Starting Haydock’s recursion from the initial state $\ket{\mathcal{K}_{0}}=\ket{\phi_{0}}=\ket{0,0}$ (electron at site $i=0$ having energy $\epsilon_{0}$ and no phonon modes excited), one then finds

[TABLE]

Here states in the energy independent tight-binding representation of $H_{IDMFT}$ are labeled $\ket{x,n}$ with $x$ being the lattice coordinate and $n$ being the phonon number (see Fig.4). Projecting on Eq.A with $\ket{\mathcal{K}_{0}}$, respectively $\ket{\mathcal{K}_{1}}$ one then finds $a(0)=\epsilon_{0}$, $b(0)=m$, and $\ket{\mathcal{K}_{1}}=\ket{1,0}$. In the second recursion step one then finds

[TABLE]

The set of new recursion coefficients then are $a(1)=\epsilon_{1}$, $b(1)=\sqrt{t^{2}+b^{2}_{\Sigma_{0}}(0)}$ and the new wave function reads

[TABLE]

This iterative three-term recurrence procedure is repeated until all recursion coefficients a(n), b(n) and Krylov subspace wave functions $\ket{\mathcal{K}_{N}}$ are determined.

To conclude, apart from $\ket{\mathcal{K}_{0}}$, $\ket{\mathcal{K}_{1}}$, which correspond to a electron localized at site $i=0$, respectively $i=1$ with zero phonon modes excited, all other vectors $\ket{\mathcal{K}_{N}}$ represents excitations of the many-body system (electron and phonon modes) that progressively spread away from the interface into the acceptor (see Fig.4).