Inelastic effects in electron transport studied with wave packet propagation.
Abstract
A timedependent approach is used to explore inelastic effects during electron transport through fewlevel systems. We study a tightbinding chain with one and two sites connected to vibrations. This simple but transparent model gives insight about inelastic effects, their meaning and the approximations currently used to treat them. Our timedependent approach allows us to trace back the time sequence of vibrational excitation and electronic interference, the vibrationally introduced time delay and the electronic phase shift. We explore a full range of parameters going from weak to strong electronvibration coupling, from tunneling to contact, from onevibration description to the need of including all vibrations for a correct description of inelastic effects in transport. We explore the validity of singlesite resonant models as well as its extension to more sites via molecular orbitals and the conditions under which multiorbital, multivibrational descriptions cannot be simplified. We explain the physical meaning of the spectral features in the second derivative of the electron current with respect to the bias voltage. This permits us to nuance the meaning of the energy value of dips and peaks. Finally, we show that finiteband effects lead to electron backscattering off the molecular vibrations in the regime of highconductance, although the drop in conductance at the vibrational threshold is rather due to the rapid variation of the vibronic density of states.
I Introduction
The importance of inelastic effects in electronic transport in molecular junctions is widely recognized and it is a rich, active research field. Several recent reviews on the topic give a clear idea of its breadth Galperin2008 (); Galperin2007 (). Advanced new experimental techniques show that electrons transport through a fewatom system is strongly dependent on the vibrational degrees of freedom of the system. As miniaturization decreases device sizes, the role of atomic vibrations needs to be considered in the device functionalities. This dependence has been shown to profoundly alter the device behavior: from new channels of conduction Berthe () to heating of the junction Neel (). It is then important to understand the role of inelastic effects and the parameters that control them.
A large body of research has been devoted to inducing controled inelastic effects. A recent review article Ho2003 () shows that inelastic effects in a tunneling junction can be used to chemically analyze the molecules at the junction by inelastic electron tunneling spectroscopy (IETS) Stipe (), to induce reactions by displacing atoms Lee (), and to use molecular conformational changes as a switch for possible devices Gaudioso (). In this way, the scanning tunneling microscope (STM) induces the reaction and also detects its product, a dramatic example is the modification and detection of trans2butene on Pd (110) Kim (). The effect of local vibrations during electron flow has also been revealed in careful experiments of photon analyses Hofluorescence () yielding an unprecedented insight in vibrational dynamics. Not only has the tunneling junction of an STM been used to explore the coupled electronvibration dynamics, also point contact spectroscopy has been used in molecular wires Ruitenbeek (); Agrait (), where the conductance showed drops at the molecular vibrational onset.
These many experimental results call for important theoretical development. Indeed, the last years have seen several theoretical works spanning most of the experimental systems: from tunneling inelastic spectra Mingo (); Lorente (); Ueba () to point contact spectroscopy Todorov (); Frederiksen (); Solomon (), from systematic approaches studying different parameters Galperin2004 (); Galperin2006 () to lowestorder expansion with abinitio parameters Paulsson (); Troisi (). A thorough review on methodology and results can be found in Ref. Galperin2007, .
A quantitative description of inelastic processes is mandatory in order to assess the relevance of the different ingredients characterizing electron transport on the atomic scale. Recent developments in abinitio calculations together with transport calculations permit us to grasp the essential parameters and, eventually, produce predictive calculations. Yet, the typical abinitiobased calculations are on the onehandside heuristical, because a well established dynamical theory of electron transport is yet to come Prodan (), and on the other handside, they are complex and difficult to interpret. Most abinitio approaches use ground state density functional theory with nonequilibrium Green’s functions (NEGF), this combination is not justified, and the results have the full complexity of NEGF. As signaled in Ref. Prodan, , new methodology to treat quantum transport will be probably based in time dependent density functional theory (TDDFT). Recent results show that it is possible to treat nuclear (semiclassically) and electronic (quantally) degrees of freedom within TDDFT to treat the nonadiabatic transport problem Almbladh (); DiVentra (). New developments go a step further in the treatment of correlated electronionic dynamics McEniry (). But timedependent methodology can have other benefits beyond its correctness. In particular, it can be used to develop a physical picture of the electron diffusion process, in this way yielding complementary information to the more involved NEGF approach. Timedependent approaches can also have interesting numerical behavior. Indeed, electron transport treated with the shortiterative Lanczos method SIL () has a quasilinear scaling for sparse Hamiltonians.
In this article, we explore inelastic effects in electron transport by means of electronic wave packet propagations in an idealized atomicsize system. We consider tightbinding chains connected to one and two vibrating electronic sites. These vibrating sites can hold one and two nuclear modes that are coupled to linear order in the nuclear displacement with the electronic degrees of freedom. Despite the simplicity of this model system, the main oneelectron ingredients are included, and the time resolved solution permits us to have insight different from the perturbationtheory Green’s functions results. Similar treatments have already been performed for the case of electronmolecule collisions Gauyacq () and for inelastic effects in transport Bringer ().
The calculations presented here are distant from the experimental situation because the model system is very simplified and because manyelectron problems are absent. Indeed, recent treatments show the richness of effects associated with the manyelectron aspects of the problem Hyldgaard (), as well as the nonequilibrium manyphonon problem Flensberg (); Mitra (). Despite, the absence of these very interesting ingredients of inelastic transport, our calculations can help in understanding inelastic effects because there are situations in which oneelectron transport is justified even in the presence of inelastic effects Imry (). To a certain degree, our calculations are equivalent and complementary to those of H. Ness Ness2006 () the main difference being that a timedependent approach is adopted in the present study.
Ii Timedependent wave packet propagation
Stationary electron transport does not need a timedependent description. However, insight on vibrational excitation processes during the current flow can be gained by timedependent calculations. Numerically, a timedependent description can benefit from the quasilinear scaling using sparse Hamiltonians. As in stationary descriptions, the bottleneck of the calculation lies in the matrix time vector product of the Hamiltonian acting on a system’s vector. Hence, a timedependent approach can yield especial insight in the actual transport process in a realistic nanoscale system by using a localized basis set that leads to a sparse Hamiltonian. This strategy seems to be particularly appealing for the description of conduction in nanowires Troels () and nanotubes Stephan (). In this section, we explore the general timedependent methodology using wave packets, beyond the Kubo linear theory of Refs. Troels, ; Stephan, .
A particularly interesting timedependent approach to have information on the full electronic trajectory is the shortiterative Lanczos (SIL) wave packet propagation SIL (). The initialvalue problem is solved by applying the evolution operator to the initial wave function, so that the solution reads, , where is the Hamiltonian of the system under consideration (atomic units are used throughout, , unless otherwise specified). The use of this equation is inconvenient in the case of large matrix dimensions because it implies a diagonalisation. Instead, we prefer to use a numerical approximation which consists in successive infinitesimal evolutions of the wave function, . In order to have an efficient implementation, the SIL method truncates the Hilbert space to a subspace spanned by a few vectors. The computation of this truncated subspace is the bottleneck of the calculation because the matrixtimevector product to generate the subspace is performed on the total dimension of the problem. Despite the truncation, the richness of the full Hamiltonian spectrum is recovered by repeating this operation for all the different time steps and propagating in this way the initial wave packet.
Let be a starting vector, the Lanczos propagation method consists in the construction of a Lanczos matrix following the recursion relation Lanczos (); Cullum (),
(1) 
where are the new diagonal matrix elements of the new tridiagonal matrix and are the new upper and lower diagonal matrix elements. Since the recurrence relation is started by , initial wave function, , then and . The vectors are called Krylov vectors, and define the Krylov space of order , by spanning the subspace given by , namely
(2) 
The Krylov vectors are orthogonal by construction, they may be normalized, and define thus a basis set in which the Hamiltonian can be expressed. In this basis, the Hamiltonian is tridiagonal. Additionally, the order can be much lower than the initial matrix dimensions. Consequently, the Lanczos matrix can be easily diagonalized using conventional algorithms. Regarding the Krylov vectors, it is known that for relatively large values of , their orthogonality can be lost, this is the reason why we have set sufficiently small time steps to keep a rather small order (). Its value can be changed at each step in order to optimize the performances of the calculation, lowering the computational cost as decreases. An excellent description of this dynamical control of the accuracy of the Lanczos propagation method can be found in Ref. Andrei, .
Another important feature is the undesirable effects provoked by the finite size of the propagation grid. To eliminate artificial reflections of the wave packet at the boundaries, we use a parabolic absorbing potential. In order to account for the errors which may be introduced by the reflections at the boundaries in the presence of the absorbing potential, we calculated the reflection coefficient with a broad wave packet. We found that less than 0.01% of the wave packet was reflected at any energy.
To compute the energyresolved transmission or reflection coefficients we use the virtual detector technique Detectors () which consists in the evaluation of the wave function some sites after the region where the electronvibration interaction takes place. By a timetoenergy Fourier transform, transmission or reflection coefficients are obtained. In the inelastic case, this should be done for each vibrational state of the total wave function. Hence, for a 1D tightbinding chain, the partial transmission, , is given by
(3) 
where we have assumed zero temperature and an initial wave packet in the vibrational ground state, . is the energyresolved wave function in the vibrational state , after the interacting region. The detector has been located on site . is the wave function computed with the same virtual detector but considering a bare Hamiltonian containing no vibrational degrees of freedom nor elastic defects.
To analyze quantities such as total transmissions, especially for nontrivial Hamiltonians which may have two or more vibrating electronic states, we have computed the density of states projected on any state, for instance, on any linear combination of the tightbinding states. Consider a state , with a particular electronic state and its vibrational state, evolving with Hamiltonian . The density of states projected on reads Tannor (),
(4) 
Iii The physical model
The FröhlichHolstein model Holstein () is used to represent the combined electronvibration system: the electronvibration coupling is assumed to be linear in the normalmode coordinates. If we assume only one mode of vibration, and a single site, the Hamiltonian reads,
(5)  
where and are the fermionic operators which create and anihilate an electron in state . Similarly, and are the bosonic operators that create and anihilate a quantum of energy in the considered vibrational mode. The first term in the Hamiltonian (III) refers to the site where the interaction takes place, it has an onsite energy of . The second describes the energy of site of chain . For simplicity we will be dealing with only two chains (: left and right). The third term describes the couplings among sites of chain . Here, we just consider nearest neighbours. The fourth term is the coupling between the two semiinfinite chains and the state of energy . The fifth term is the energy of the harmonic oscillator. Finally, the last term describes the electronphonon interaction. In this singlestate singlemode model, is only a scalar which represents the strength of the interaction. It is called the electronphonon coupling. In the present work, we have explored both an electronic single site coupled to vibrations, and a double site. In the latter case, will be a matrix, as presented in section V.
Using a tensorial product description of the electronic and nuclear coordinates, the full Hamiltonian can be expressed in matrix form:
(6) 
where we have used a block representation in the vibrational basis . The diagonal elements are electronic Hamiltonians, which define the propagation of a wavepacket in a vibrational subspace. For clarity we give their tightbinding representation,
(7) 
where labels a particular vibrational state in the harmonic approximation. The are offdiagonal matrix elements which connect the nearest neighbours sites inside the left and right chain. The end play the same role as , they connect the chains to the site which has onsite energy . The diagonal term accounts for the energy the electron must exchange with the vibrational degrees of freedom of the system. If it propagates from one vibrational state to another, it must lose or gain , the energy quantum of the vibration.
The matrices in (6) couple the Hamiltonians in the different vibrational subspaces. In the case of a singlesite impurity, is essentially a sparse matrix with the same dimensions as , where only one single element is nonzero, the one connecting the diagonal matrix elements of energy and . In Hamiltonian (6), the factors that multiply come from the matrix representation of the bosonic operators, they correspond to the factors which appear in the well known relations and . We note that Hamiltonian (6) is truncated, the number of vibrations that are considered in the calculation is in this example. This number may be sufficiently large to represent suitably the vibrational space. A detailed study of the convergency of the calculation is presented hereafter.
Let us consider a molecule between two electrodes, modeled by a single state connected to two chains. In this singlesite case, Y. Meir and N. S. Wingreen Meir () have shown that the current can be expressed as follows,
(8) 
where is defined as a function of the couplings to the right and left leads, . is the retarded green function of the molecule, and is the Fermi distribution of the left (right) lead. This formula applies in the case where the couplings to the leads only differ by a constant factor, , such that . This is always so in the wideband limit for the singlesite case. This is not the case beyond the wideband limit.
The equation above can be viewed as the integral of a quantity that we identify as a transmission, multiplied by an energy window given by the difference of the Fermi distributions of the leads, meaning that a current will flow if both the transmission is nonzero and the voltage applied between the electrodes is sufficiently large. Following the derivation by H. Ness Ness2006 () at the zero temperature limit, the transmission reads,
(9) 
where is the projection of the Green function in the vibrational subspace, . This means that the transmission is related to the projected density of states, Eq. (4), in such a way that, if we consider the wideband approximation, where the coupling is independent of energy, the transmission is proportional to the density of states projected in the subspace. By using the optical theorem one can explicitly retrieve the vibrationally excited states in the inelastic current Ness2006 ().
In the case of the wave packet calculation, we can extend these equations for the calculation of the electronic current. We approximate the conductance assuming that the Fermi level only enters to define possible final electronic states (otherwise the calculation is fully oneelectron, an in general we will not consider a Fermi level)
(10) 
where the factor is the conductance quantum in atomic units. The terms between brackets are introduced to account for the opening of vibrational channels since they contain the Fermi distribution functions of the left electrode, and the right one, , at zero bias voltage in the present case. Actually, when electrons do not have sufficient energy to deposit it into the vibrational degrees of freedom of the system, only the is to be considered in the conductance. The same holds if the energy of the incident electron is sufficient to excite one vibration, in this case we consider the sum of and . Nevertheless, the calculation should be performed with a sufficiently high value of in order to take into account the influence of closed channels in the conductance. Finally, the current as a function of the voltage, , can be written as the integral over energy of the conductance,
(11) 
The terms between brackets ensures that the transmitted electrons go from occupied states to empty ones where the voltage dependence is included.
The treatment presented here is complementary of the one by H. Ness Ness2006 (). Indeed, the physical model is the same one, the difference is the solution method. As in the case of Ref. Ness2006, the present approach is singleelectron, neglecting both electron occupation and electronelectron correlation effects.
Iv Singlesite resonant model results
Let us assume a single state coupled to the two 1D tightbinding chains of the above model. Figure 1 shows the results of a wave packet propagation. Figure 1 (a) is the elastic wave packet, , because we assume that no vibration is initially excited in the system and the temperature is zero. The present system is two 1D electodes symmetrically coupled to the vibrating site. Hence, in Fig. 1 (b) we see that only the reflected wave packet is different from the transmitted one in the elastic channel, , while it is identical for the inelastic channels: the electron reaches the active site and populates the ladder of excited vibrational states. In each state the electron has a finite and identical left and right transmission.
A simpleminded picture of the structure appearing in the transmission as a function of incident electron energy could be that there is a series of N resonances displaced by the phonon energy, in agreement with the scheme of Fig. 2. At a given inital electron energy, the wave packet probes the resonant electron site if the energy is within of the resonant site energy. The results is a series of equidistant peaks spaced by . This is plotted in Fig. 3 (a). There, the transmission for the single electronic state without electronvibration coupling is depicted in a dotted line, and the full line corresponds to the case when the electronvibration coupling is included.
However, the physical picture is more complex than just a series of equidistant resonances. In order to understand the appearance of the vibrational sidebands we divide the transmission according to the final state of the vibrator, Fig. 3 (b). This yields information on the vibrational state once the wave packet has propagated through the resonant site. The result is that each transmission curve for a welldefined final vibrator state displays a similarly rich peak structure. The above picture, Fig. 2, has to be changed: there are complex vibrational pathways, in which different vibrational states are probed before the system is left in a singly welldefined vibrational state. Given the coherence of the electron propagation, the different pathways can interfere and will give rise to Fano lineshapes in the transmission function Durand ().
iv.1 Vibrational state population sequence: electron coherence
Wavepacket dynamics can probe the different structures appearing in the transmission function in order to yield information on the actual interference process. Let us take a vibrational ground state, broad electronic wave packet energetically centered about the first maximum of the full transmission function, Fig. 3. A spatially broad wave packet is almost a planewave and synonymous of a monochromatic packet, hence we can be sure to probe only the structure of the first maximum.
Figure 4 shows the modulus square of the electronic wave function for the first three vibrational levels on the vibrating site as a function of time. We first notice that the population of the vibrations are sequential: the peaks are slightly shifted in time as increases. The linear electronvibration coupling forces this sequential population since is changed in steps of one quantum of vibration, Eq. (III). The populations show a time modulation, in the present case, we see that the modulation frequency of the wave packet is roughly twice the modulation of the ground state and of . This is easily explained by considering populating and depopulating the level by populating the and the . After a certain time the level can depopulate again in the as well as the . Despite the simplicity of the electronvibration coupling and the electronic model, the final population dynamics depend on the wave packet energy, the strength of the vibrational coupling and the electron lifetime at the site, leading to nontrivial dynamics. The lifetime of the site is fixed by the hopping matrix elements and , Eq. (7).
The Fano profile characterizing the transmission functions, Fig. 3, can then be explained by the interference of several of these vibrational pathways. Indeed, asymmetric Fano profiles are the rule in Fig. 3 (b). For the case of the total transmission, the addition over final vibrational states smear out the different profiles, rendering more difficult the determination of an asymmetric Fano profile in the peak sequence.
iv.2 Time delay and phase shifts
Spatially narrow wave packets contain most of the energy components needed to span the spectral region of interest, i. e. the main peaks of Fig. 3. The Fourier coefficients in the energy domain can be evaluated to obtain the phase of each component and hence the phase shift between Fourier components. The derivative of the phase shift with respect to energy yields the time delay of the Fourier component in the wave packet Gauyacq (). Hence, we can analyze the effect of the vibrational excitation of the electronic site by means of the time delay imposed on the propagating wave packet.
Figure 5(a) shows the phase shift for the elastic component of the wave packet (the component). The full phase shift over all of the vibronic peaks is , this is the phase shift of the electronic resonance. The phase shift of each individual contribution to the transmission is more complex. We see that the phase shift varies rapidly after each maximum, very much indicating a suite of BreitWignerlike resonances. Indeed, the the phase shift for a BreitWigner resonance of FWHM centered at is given by Gauyacq ():
The time delay, , is given by the energy derivative Gauyacq () of the phase :
In the case of a single BreitWigner resonance, the time delay is then
On resonance the time delay is just , yielding direct information on the resonance width, . In the present case, the time delay, Fig. 5(b), gives some definite interpretation: the time delay is maximum once the wave packet has encompassed one of the vibronic resonances. Hence, just above the resonance, the electron is deterred in its propagation by interaction with the vibration. The vibration contributes to the partial width of the electronic resonance, however the total width is independent of the vibration Wingreen (). A spatially narrow wave packet, will have the phase shift and time delay of an electronic resonance regardless of the existence or not of the vibrations. We also find negative timedelays, that correspond to drops in the phase shift. Electrons can hence be expelled from the resonance more easily in the presence of vibrations at certain incident energies. This behavior is due to the interference between two consecutive vibronic resonances. As emphasized in Ref. Gauyacq, , the time delay has to be interpreted with care in the case of a multichannel problem such as the present one. The time dependence is indeed complex and the time delay is just a number that cannot summarized the full intereference pattern.
The interferences among vibrational paths are ubiquituous in all present results. Their effects can be seen in the Fanolike lineshapes of the tranmission function Fig. 3, in the oscillating population of vibrational states with time, Fig. 4 and in the rapid changes of the phase shifts and in the negative time delays, Fig. 5.
iv.3 Time scales
To develop an understanding of the vibrational excitation process, many authors resort to comparing the different timescales in play Galperin2008 (); Galperin2007 (); Baratoff (). This is a very appealing approach because it permits us to rationalize the excitation process and the excitation regime of different systems.
Let us define as the lifetime of an electron in the molecule, which is given by the inverse of the resonance width, . Typically, for chemisorbed molecules, eV, which means electron lifetimes in the subfemtosecond range. In order to estimate the change in conductance due to an inelastic process, we can compute the inelastic fraction of electrons passing by the molecule. This can be estimated by computing the ratio of the excitation time to the first quantum of vibration, , and the lifetime of the electron in the molecule, . Let us estimate the magnitud of the electronvibration coupling, , to attain a measurable change in conductance (larger than 1%). This is the regime where perturbation theory is valid. In this case, the excitation process has been classified as sudden in the electronmolecule collisions literature GauyacqLibro ().
In a sudden process, the molecule is assumed to be very briefly in its negative ion potential energy surface (PES). This can be indeed very brief as has been said above. As in Ref. Berthe, , let us assume that the PES is basically parabolic, but displaced with respect to the neutral PES Negative ():
(12) 
where is the spring constant of the PES related to the mode frequency by where is some reduced mass (this is easily generalized to the case of multinuclear modes). Here, is the displaced center of the negative PES. When the electron flows through the molecule, the molecule is suddenly in its negative PES. Hence, the nuclei experience a force given by
Then, they acquire a speed, , of the order of
and . Hence,
where is the electronvibration linear coupling, of Eq. (III) with coming from the second quantization of the displacement . The speed gained by the nuclei at an excitation of one quantum of vibration near is . In this way, we find an upper limit for :
(13) 
This simple estimation shows that for strong electronvibration coupling, vibrational excitation is unavoidable, where strong means larger than the moleculeelectrode coupling.
Weak coupling is then the regime when . We can use Fermi’s golden rule to estimate the vibrational excitation rate Mats ():
Assuming typical IETS branching ratios , leads to
therefore the coupling becomes .
For chemisorbed molecules this is in the range of 0.01 eV.
This coupling is easily found
in a large class of molecules, and therefore IETS is a feasible spectroscopy. This
is a simple estimation showing that rather small ( eV)
can have a measurable change in conductance for molecules adsorbed on metallic leads.
One more time scale is fixed by the vibrational frequency, . As discussed in Ref. Gauyacq, , there are three basic regimes that can be distinguished:

the negative ion is long lived and the molecule can vibrate in the new state,

the negative ion is short lived and we are in the above sudden regime,

strong interference effects appear due to the nuclear motion.
The time scale given by , defines then the type of electronvibration interaction that will result. In the first case, the electron has ample time to interact with all the nuclear degrees of freedom and depending on the electronvibration coupling, , vibronic signatures can appear. This has been shown in the case of STM studies of electronic states on surfaces Berthe (); Repp () with a band gap leading to small . The second case corresponds to the IETS case discussed above. The third case has been studied in the gas phase (see for example Ref. Gauyacq, ) but we are not aware of any report on the consequent boomerang effect in the electron transport regime.
iv.4 Finiteband effects
The wideband approximation simplifies the expression for the transmission that becomes analytical Wingreen (). One of the results of the wideband approximation is that the transmission function and the spectral function or PDOS, Eq. (4), have the same behavior with the electron energy. This is definetly not the case when finiteband effects are considered, indeed, the energy dependence of the coupling to the electrodes leads the transmission function to have a different energy behavior than the site spectral function, Eq. (9). When the initial level, is near the center of the band, the conditions for the wideband limit are easily attained. As the level approaches one of the band edges, becomes strongly dependent of energy as well as the center of the transmission distribution, Fig. 3.
One of the apparent features is that the transmission peaks become thinner and better defined. As approaches the band edge, many inelastic channels start closing because the final electron energy falls out of the electron band, leading to a substantial decrease in the transmission peaks. This is clearly seen in Fig. 6. The parameters of the transmission calculated in Fig. 6 are exactly the same as for Fig. 3 (a) except that is 0.95 eV instead of 0.7 eV. The bottom of the band is at 1 eV and the quantum of vibration is eV. The singular aspect of the peaks is enhanced by the 1D character of the electrodes. Indeed, when the final energy of the electron approaches the bottom of the band a van Hove singularity appears in the density of states making the transmission more singular. In more realistic cases with 3D electrodes we expect less singular transmission peaks.
For an electron resonance approaching the top of the band, the situation is rather different. Here, the higher vibronic side bands dissapear but the convergency remains the same because no channel is closed, leading to lowenergy side bands of the same width as for the case of the resonance at the center of the band.
iv.5 The dynamical polaron shift
The polaron shift is defined as the energy displacement of the first peak in Fig. 3 (a) with respect to the dashedline peak. The polaron shift is related to the appearance of a vibronic structure. As we showed above, in the timescale discussion, 2 conditions need be satisfied. Namely, must be sizeable and . This leads to considering Hyldgaard () the parameter . When the polaron shift becomes measurable.
Hyldgaard et al. Hyldgaard () show that the polaron shift critically depends on the initial occupation of the vibrating electron site. When the electron site is occupied, the transmission function shows a main peak displaced by . At halffilling the displacement is and for an empty site the polaron shift is .
We can give a direct interpretation of the polaron shift by using the above parabolic model, Eq. (12). Following the above discussion, the polaron shift for the empty site is which is equal to . The neutral PES is given by , hence the polaron shift is exactly the energy difference between the neutral PES and the negative one at the nuclear coordinate at which electron capture takes place. Hence, the polaron shift is the energy gain in the formation of the negative intermediate because the electronvibration coupling permits the electron to probe the negative ion PES. In the absence of coupling, the nuclear coordinates do not evolve, and the negative ion resonance is just a simple lorentzian.
In order to probe the polaron peak, the nuclear wavefunctions in the two ground states (neutral and negative) need to overlap. The overlap is large at weak electronvibration coupling, hence a welldefined peak appears, just shifted by the energy gain. As the coupling increases, the overlap diminishes, leading to a decrease of the polaron peak. In the limit of strong coupling, the polaron peak is basically zero and a rich structure of evenly distributed peaks of the negativeion vibrational states is apparent, with a gaussian envelop Repp (); Mahan ().
The polaron shift is a good test for the convergency of numerical calculations with the number of vibrational states included in the expansion of Eq. (6). Figure 7 shows the convergency of the polaron shift with the number N of phonons included in the calculation. A correct value is obtained () at about N=5 phonons. Following the results of the calculation of Hyldgaard et al. for an empty state, the polaron shift would be with our parameters ()
We also show in Fig. 7 the interpeak distance as a relevant quantity regarding the convergency of the calculation. We observe that the peaks of lower energy tend to converge with less phonons than the peaks of higher energy, which need a higher number of phonons to be placed at the one with respect to the other. Trivially, higherorder peaks need higher , otherwise the peak may not even appear in the calculation.
The polaron shift is also very sensitive to the different approximations that are used when evaluating inelastic effects in electron transport. Indeed, the selfconsistent Born approximation (SCBA) Hyldgaard () yields wrong polaron shifts. Such a study is perfomed in Ness2006, . There, it is shown that the SCBA is equivalent to neglecting the factors of Eq. (6), this leads to wrong interpeak distances. Nevertheless, the SCBA captures much of the physics of inelastic processes and is an allorder theory, becoming reliable enough and very interesting for the evaluation of inelastic processes in a wide range of problems. In the particular case of weak electronvibration coupling, and the SCBA is virtually exact.
V Twosite resonant model results
Most of the literature devoted to transport in the presence of electronvibration coupling is based on the singlesite model. However, Bringer et al. Bringer () showed that the singlesite case has a behavior that cannot be extended to more realistic systems in which several electronic states are coupled with vibrations. Here, we will analyze when one can reduce the problem to the singlesite case and when not. In the same way, we will show that when several vibrations are involved, the neglection of some of the vibrations can lead to qualitatively wrong results.
v.1 Twosites and one vibration
Our previous onesite model is extended to have two electronic sites connected to a vibration. This model has been recently explored in the context of abinitio based calculations of inelastic transport between two pyramidal electrodes Frederiksen2007 (). The simplified twosite model permits us to understand the more complex abinitio results. The model is given again by Hamiltonian (III) and (6) where and are matrices.
Here, one only vibration is assumed to interact with the electron flow. In the case of weak coupling Frederiksen2007 (), this assumption is justified because Eq. (III) is truncated to the first phonon of each mode and the inelastic effects of each mode are separable. In the next section, we will see that when convergency in the number of phonons needs more than one excited vibration, all of the vibrations need to be considered at once.
In this onedimensional model, there are two only possible modes Frederiksen2007 (). We will call them the centerofmass mode (CM), and the antibondlength mode (ABL). The CM mode means that both sites displace in phase, hence the electronvibration coupling matrix have nonzero onsite elements (the diagonal):
(14) 
and the second element is negative in order to account for the sign of the displacement of each site. In the ABL mode, the sites are moving out of phase, corresponding to an internal stretch mode. Hence, the electronvibration matrix becomes
(15) 
Figure 8 shows the transmission for the ABL and CM modes, independently computed and together (to be analyzed in the next section). These figures show series of vibronic peaks difficult to analyze.
In order to analyze the electron transmission, Fig. 8, it is more convenient to use projected density of states (PDOS) of the full vibronic structure on molecular orbitals, Fig. 9(a) for the ABL mode and (b) CM mode. These orbitals are linear combinations of the two sites:
and
that diagonalize the uncoupled twosite Hamiltonian with eigenvalues and , respectively, where is the level energy for each site and is the hopping matrix element between sites.
In this new basis set, the above coupling matrices, Eq. (14) and (15), become:
(16) 
and
(17) 
that hint at the interpretation of the above vibronic peaks. In the case of the CM mode, the and orbitals are coupled via the vibration, while in the ABL case, the molecular orbitals are not mixed by the electronvibration interaction. This case can then be interpreted as two singlesites connected to a vibration, and hence two vibronic sequences are associated with the spectral feature of at and the spectral feature of at . Since, the effective electron vibration coupling is for , the number of peaks and the general vibronic sequence corresponds to stronger coupling than the vibronic sequence of the peak. This is clearly seen in the PDOS on and , Fig. 9, where we see that the vibronic sequences with and character are energetically localized near the original molecularorbital peaks.
The CM mode mixes and , hence we obtain vibronic peaks that are shared in the PDOS over both molecular orbitals.
For small enough couplings, we can estimate the vibronic structure keeping just in the vibrational part. Hence we can diagonalize Hamiltonian (III) for the sites connected to the vibrations. This can also be expressed by a hybridization scheme, Fig. (10). We realize that the new peaks correspond onetoone to the peaks found in the PDOS, and that they have contributions in different ratio from the and orbitals according to the matrix elements of M. Depending on the magnitude of more or less phonons, , will be needed to converge the vibronic sequence Mahan (), and more or less elements will be included in scheme (10). In order to obtain the type of diagram (a) or (b), only knowledge of the symmetry of the system is required.
In the present case, while the CM mode mixes and orbitals efficiently, the ABL mode does not mix them. Hence, only this last case can be understood by a singlelevel model.
It is interesting to study the time dependence of the electron occupation of both sites depending on the considered mode. Given the “shuttlelike” motion of the CM mode, one could expect that both sites populate at the same time when they approach together the left electrode and then depopulate when they approach the right electrode. On the other hand, the ABL motion is a stretch mode and hence, one site would populate while the other one depopulates.
However, in Fig. 11 (a) we see that the ABL mode leads to population of both sites at the same time, while in presence of the CM mode, Fig. 11 (b), the population sequence is shifted and dependent on time. In this case, the wave packet has been centered about the energy and its energy span is smaller than the difference between and . The reason for this behavior is that one should think in terms of molecular orbitals rather than site orbitals. In this case, the ABL motion does not couple the and orbitals, hence, the electron populates one of the molecular orbitals, , at one time and interacts with the vibration in the same way as in the case of the single site. Instead of a single site, we have a single level, otherwise there is no physical difference. Let us approximate the wave packet by the evolution on the two sites:
without any time dependence (except the one due to the coupling to the electrodes that we have neglected in this simplified discussion) because is basically an eigenstate of the Hamiltonian. This explains Fig. 11 (a).
In the case of the CM mode, we saw above that a singlesite analogy cannot be applied and we have the two molecular orbitals involved in the wave packet propagation. Let us again approximate the wave packet by the evolution on the two sites:
where a clear time dependence subsists. In good agreement with Fig. 11 (b).
v.2 Twosites and two vibrations
When more than one phonon per mode is needed to be converged with respect to the electronvibration coupling, the inelastic transmission cannot be separated in additive contributions from different modes and all of them have to be considered at the same time.
In the above case, the truncated Hamiltonian reads now:
Here the subindex refers to the ABL mode and to the CM one. At intermediate coupling, we need some 3 phonons to converge the transmission. Hence, the Hamiltonian matrix in this basis set is composed of blocks, each block contains the infinite number of electron sites, that we deal with as above.
The inclusion of several modes changes the transmission from the superposition of peaks coming from each mode because the spectral weights change and the peaks shift. This is seen in Fig. 8. There the peaks of the full transmission are shifted with respect to the sum of peaks from the CM and ABL modes. This can be understood in terms of the different Hilbert space that is considered as new modes are included. Indeed, since the convergency in phonons, , means the shift of peaks, the transmission also needs to be converged with respect to modes. In terms of Green’s function this is easily understandable because there are terms in the selfenergy stemming from all modes.
We can improve our previous hybridization diagram to describe the peak structure, by including all modes. This type of correlation diagram has been termed “progression of progressions” and used in the literature to explain the vibronic sequence of C and naphtalocyanine molecules Ho2005 (); Ho2007 (). One word of caution is important: the knowledge of the actual coupling matrices is needed in order to perform the correct hybridization of orbitals. The hybridization scheme can be obtained by symmetry arguments alone, however the separation and strength of the vibronic peaks need a quantitative evaluation of the matrix elements, which is increasingly difficult with the number of involved modes.
Vi The meaning of dips and peaks in the
The increase or decrease of conductance over the vibrational threshold has been used to determine the occurrence of vibrational excitation during electron flow through molecular electronic states Stipe (); Agrait ().
A simple explanation for the increase of conductance in the tunneling regime was already advanced in IETS of metalinsulatormetal junctions Hansma (). The conductance increases because new conduction channels become available above a vibrational threshold. This interpretation is correct when the vibrational side bands of the electron transmission lie beyond of the electrode’s Fermi energy. As we have already seen, opening a new channel means that we have to add the contribution to the contribution of the transmission. Hence, the will present a positive peak at positive voltage, and negative at negative voltage Hansma ().
However, we can also have decreases of conductance. This is experimentally found in tunneling Ho_O2 () and in contact regimes Agrait (). In the present theory, one always adds a positive contribution to the transmission above the vibrational threshold, but the change in conductance is given by the slope of the transmission, and this can have a rapid variation in the presence of a sideband. In the case where the sideband is in the tail of the main peak, and the width of the main peak is in the order of the vibrational frequency, then the slope will change from a smoothly varying one to the faster varying sideband slope, giving rise to a negative at positive voltage. Let us give two examples of the above cases.
vi.1 Tunneling regime
In the case of the excitation of the CO frustrated rotation mode Lauhon (), the and the molecular orbitals are coupled via the electronvibration coupling. This is easily seen by symmetry arguments Lorente2005 (), since the mode is antisymmetric with respect to the planes containing the molecular axis and the nonzero matrix elements will couple a symmetric orbital () with an antisymmetric one (the ) Lorente2000b (). Assuming weak coupling, we can neglect the effect of all other modes and estimate the change in conductance by using the twosite model coupled to a vibration that we presented above. Figure 12 shows the result of the conductance, Eq. (10), and in the corresponding . We see that the Fermi energy plus for the frustrated rotation, gives the threshold where the transmission for enters, giving a clear discontinuity in the conductance (positive contribution) and hence a positive peak in the . This is the case in metalinsulatormetal junctions Hansma (). The twolevel model leads to different inelastic efficiencies at positive and negative voltage bias.
However, the case of IETS of O on Ag (110) Ho_O2 () is different. Experimentally, dips are found instead of positive peaks in the . Let us consider the excitation of the internal stretch mode. By computing the PDOS on molecular orbitals Olsson (), we know that the perpendicular to the surface is close to the Fermi energy, and that tunneling takes place through it Olsson (). The mode has the same symmetry as the molecular orbital and the diagonal matrix element will be different from zero. We are hence in a case that can be approximated by the singlesite model. Figure 13 shows the result. We have located the at the Fermi energy as the PDOS on molecular orbitals Olsson () seems to indicate. The orbital width is taken as 0.1 eV and the frequency energy is 0.08 eV. We see that the transmission function is an asymmetric peak because the sideband is inserted in the tail of the main peak. We also include the contribution of above the Fermi energy plus the frequency in order to calculate the , as for the CO case while we prefer to plot the complete oneelectron transmission because in the present case the initial level occupation matters and we cannot treat it in this oneelectron approach. Instead, we realize that the see the change of slope due to the side bands in the and contributions. The gives a dip corresponding to the change of slope of the contribution. Hence, the elastic component of the transmission contains the information of the IETS. The decrease of conductance in this case is due to the rapid change of transmission already in the elastic channel. It is then a decrease due to the variation of the vibronic density of states.
From these models we conclude that negative peaks in the are due to sidebands in the elastic transmission. This statement is equivalent to the one found in the pioneering work by Davis Davis (). Davis realized that a decrease of conductance could be found in some especial circumstances. Namely, that the singlesite resonance was at the Fermi energy and that the resonance width was of the order of the vibrational frequency. These are the same conditions as we find. The interpretation by Davis is that virtual phonons were emitted and absorbed giving rise to an interference pattern. This interpretation is based on perturbation theory, and it just accounts for the vibration’s effect on the elastic channel. In our terms, it is just the vibrational sideband of the elastic channel. The only cases in which a vibrational sideband can yield a decrease in conductance is when the main peak of the transmission is at the Fermi energy Persson () (this is halffilling, we discuss below that we cannot have halffilling in our models), and the electronic widths are of the order than the vibrational frequency, so as to enlargened and distort the main peak by the sideband. In the case of halffilling, the sidebands are symmetric with respect to the Fermi energy Hyldgaard () giving rise to a peak at negative voltage and to a dip at positive one in the . The voltage of the dip does not exactly correspond to the vibrational frequency, since it correspond to the largest slope of the vibrational side band, not its maximum. Hence, these results suggest that the voltages at the dips are not direct measurements of the vibrational frequencies, while in the case of the voltages corresponding to peaks in the , are direct readings of the vibrational frequencies.
There is certain confusion in the literature because perturbation theory is currently used. In this case the current in the absence of electronphonon coupling is confused with the elastic current. As can be seen in Fig. 3 (a) in the absence of electronphonon coupling one gets a single lorentzian peak (dashed lines), however the elastic contribution is the curve of Fig. 3 (b). Hence, it is wrong to identify the elastic current with the one without electronvibration coupling.
vi.2 Contact regime
Paulsson and coworkers Magnus () have shown that in the case of symmetric contact to the electrodes one can continuously pass from a peak to a dip in the by increasing the transmission probability. At transmission 1/2 the threshold between both behaviors is found. In the case of contact, a large density of states is pinned at the Fermi energy. Hence, it is the parameter that controls the coupling to the leads and, therefore, the passage from peaks to dips in the change of conductance at the vibration threshold. Recent experimental evidence has been reported in Tal, .
Figure 14 shows the behavior of the electronic wave packet after collision with a 4site chain in the contact regime. The parameters are those from Ref. Thomas2004, . The chosen mode is the ABL, that corresponds to the one giving the largest change in conductance in gold monoatomic chains Agrait (). The wave packet is mainly transmitted in the elastic channel as is to be expected from the contact regime, however the vibrationally excited wave packets are mainly reflected in agreement with the interpretation of Ref. Agrait, . However, the behavior is more complex: electronic wave packets in channels with odd are reflected, even are transmitted. In this particular calculation, the reflection of the component is due to finiteband effects. As we discussed above, a finite band leads to the closing of transmission channels when inelastic effets are operative, increasing the electronic wave packet reflection. As the number of vibrational excitations increases, the closing of channels corresponds to each excitation process, implying an increase of the reflection in each case. As a consequence, even excitations mean even number of reflections, leading to an increase in transmission. Indeed, we have found this behavior for large transmission even in the singlesite case.
Some of the calculations of Ref. Magnus, have been performed in the wideband approximation, and hence the closing of channels due to finiteband effects is not present. Hence, we cannot conclude that the decrease in conductance leads to an increase in electron reflection such as the one found in this section. Instead, we think their calculations correspond to the regime where rapid changes in the vibronic density of states are found which correspond to the precedent section.
Vii Manyelectron and electronelectron efects in the presence of vibrations
Our calculations are exact regarding the vibrational dynamics during electron transport. The price to pay is the neglect of multielectron effects. The description of multielectron effects together with vibrational excitation necessarily implies approximations. A recent treatment of both effects is the one by Galperin et al. Galperin2006 (). Despite approximations, their treatment is correctly evaluating the vibrational dynamics. Indeed, we can reproduce their results, and the main difference between both calculations is an absorption shoulder appearing below the first peak in the multielectron case. According to Ref. Galperin2006, this is due to the propagation of holes coexisting with electron transport. This is absent in our oneelectron treatment.
Even in the absence of a direct electronelectron interaction in the FröhlichHolstein model used here, Eq. (III), the electronvibration coupling can lead to effective electronelectron interactions, see for example Ref. Hyldgaard, . This leads to electronelectron correlations and has extensively been studied in the literature in the socalled bipolaron problem Zhang (). This is absent in our oneelectron approach, but could in principle be treated by using twoelectron propagations.
Perhaps, more important for the description of transport in metallic systems is the absence of the initial electron occupation in the present approach. This is easily solved by usual nonequilibrium Green’s functions approaches Hyldgaard () at the price of simplifying the treatment on the phononic description. Indeed, while the present approach treats the nuclearcoordinate description exactly within the oneelectron problem and the harmonic approximation, NEGF approaches need to rely on approximate treatments Galperin2006 (); Hyldgaard (); Mitra (). This situation can be solved by using a timedependent HartreeFock approach to generate a manybody wave packet, extending the present treatment to a multielectron case. Some encouraging results in timedependent HartreeFock and beyond techniques have been presented in Refs. Nest, ; Krause, .
Viii Summary and conclusions
Timedependent wave packet propagations for the study of inelastic effects in electron transport can be very interesting because of the physical picture they permit us to develop as well as the good sizescaling properties that they can have. Indeed, in the case of sparse Hamiltonians, the scaling is basically linear with the system size, when Lanczoslike propagation methods are used.
In the present work we have analyzed the electronvibration problem when a single electronic site is connected to two electrodes. We have used a timedependent description to understand the vibrational excitation sequence, the electronic phase shift and the interference patterns. The timedependent calculation has given us access to the same quantities we can calculate using an energyresolved theory, but with the timedependent insight.
We have also studied the case of two sites and two vibrational modes, and we have realized that the twosite problem can be simplified to the singlesite one when molecular orbitals can be used. The electronvibration coupling will determine when molecular orbitals are meaningful. This is seen by changing the electronvibration coupling matrix to the molecularorbital basis set. Pursuing this idea, we have developed a simple hybridisation scheme that permits us to understand the vibronic structure of an electron transmission function in terms of the symmetry of the electronvibration couplings. We have shown, that for medium and strong electronvibration coupling all modes need to be considered.
We have applied our simple model to get some insight in the case of the vibrational excitation of CO on Cu(100) Lauhon () and O on Ag(110) Ho_O2 (). Due to our oneelectron treatment, the calculations are missing fundamental ingredients, but they permit us to explain the IETS signals in terms of vibrational channels. In the case of O, the calculation show that the elastic channel leads to a decrease in the . Contrary to perturbation theory, the elastic channel is not the electronic structure in the absence of electronvibration coupling, but the transmission function. Hence the elastic channel contains all of the information about the electronvibration coupling. In particular, the elastic contribution will contain vibrational satellites. These satellites or side bands have a strong energy dependence that naturally leads to dips in the . Indeed this means that the density of states that should be considered is the vibronic one, this density of states contain rapid variations with the energy that lead to decreases of . In perturbation theory, one needs to consider the effect of the vibration in the elastic channel to retrieve the effect of the vibration in the electronic structure and hence to take into account the vibronic density of states in the conductance.
We have also explored the case of contact with the two electrodes. In this case, it has been shown that the conductance should drop at the vibrational thresholds because the vibrations backscatter the impinging electrons. Our calculation show that for the first excited (and generally odd vibration channel, ) the electronic wave packet is reflected at the vibration. However for even the wave packet is transmitted. This behavior is due to the finite electronic band of the present calculation, that leads to the closing of inelastic channels increasing reflection. An even number of consecutive reflections lead to increased transmission.
New developments in timedependent HartreeFock lead us to think that a full timedependent treatment of the physics usually explored with nonequilibrium Green’s functions will be soon available permitting us to develop a different point of view on inelastic effects in electron transport on the atomic scale.
Acknowledgements.
Interesting discussions with A. Arnau, A.G. Borisov, T. Frederiksen, J.P. Gauyacq, E. Jeckelmann, T. Klamroth, M. Paulsson, P. Saalfrank, and H. Ueba are gratefully acknowledged. N.L. acknowledges support from the Spanish MEC (No. FIS200612117C0401). Computer resources of the Centre de Calcul de MidiPyrénée (CALMIP) are gratefully acknowledged.References
 (1) M. Galperin, M. A. Ratner, A. Nitzan, and A. Troisi, Science 319, 1056 (2008).
 (2) M. Galperin, M. A. Ratner, A. Nitzan, J. Phys.: Cond. Matter 19, 103201 (2007).
 (3) M. Berthe, A. Urbieta, L. Perdigao, B. Grandidier, D. Deresmes, C. Delerue, D. Stievenard, R. Rurali, N. Lorente, L. Magaud, P. Ordejon, Phys. Rev. Lett. 97, 206801 (2006).
 (4) N. Néel, J. Kröger, L. Limot, T. Frederiksen, M. Brandbyge, and R. Berndt, Phys. Rev. Lett. 98, 065502 (2007).
 (5) W. Ho, J. Chem. Phys. 117, 11033 (2002).
 (6) B. C. Stipe, M. A. Rezaei, W. Ho, Science 279, 1907 (1998).
 (7) H. J. Lee and W. Ho, Science 286, 1719 (1999).
 (8) J. Gaudioso, L. J. Lauhon, W. Ho, Phys. Rev. Lett. 85, 1918 (2000).
 (9) Y. Kim, T. Komeda and M. Kawai, Phys. Rev. Lett. 89, 126104 (2002).
 (10) X. H. Qiu, G.V. Nazin and W. Ho, Science 299, 542 (2003).
 (11) R. H. Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C. V. Hermet, and J. M. V. Ruitenbeek, Nature 419, 906 (2002).
 (12) N. Agrait, C. Untiedt, G. RubioBollinger, S. Vieira, Phys. Rev. Lett. 88, 216803 (2002).
 (13) N. Mingo and K. Makoshi, Surf. Sci. 438, 261 (1999); Phys. Rev. Lett. 84, 3694 (2000).
 (14) N. Lorente and M. Persson, Phys. Rev. Lett. 85, 2997 (2000).
 (15) T. Mii, S. Tikhodeev, and H. Ueba, Surf. Sci. 502, 26 (2002); Phys. Rev. B 68, 205406 (2003).
 (16) M. J. Montgomery, T. N. Todorov, J. Phys.: Cond. Matt. 15, 8781 (2003).
 (17) T. Frederiksen, M. Brandbyge, N. Lorente and A.P. Jauho, Phys. Rev. Lett. 93, 256601 (2004).
 (18) G. C. Solomon, A. Gagliardi, A. Pecchia, T. Frauenheim, A. Di Carlo, J. R. Reimers, N. S. Hush, J. Chem. Phys. 124, 094704 (2006).
 (19) M. Galperin, M. A. Ratner, A. Nitzan, J. Chem. Phys. 121, 11965 (2004).
 (20) M. Galperin, A. Nitzan, M. A. Ratner, Phys. Rev. B 73, 045314 (2006).
 (21) M. Paulsson, T. Frederiksen, M. Brandbyge, Nano Letters 6, 258 (2006).
 (22) A. Troisi, M. A. Ratner, and A. Nitzan, J. Chem. Phys. 118, 6072 (2003).
 (23) E. Prodan and R. Car, Phys. Rev. B 76, 115102 (2007).
 (24) C. Verdozzi, G. Stefanucci and C.O. Almbladh, Phys. Rev. Lett. 97, 046603 (2006).
 (25) M. Di Ventra and R. D’Agosta, Phys. Rev. Lett. 98, 226403 (2007).
 (26) E. McEniry and T. Todorov, J. Phys.: Cond. Matt. 19, 196201 (2007).
 (27) S. Mahapatra and N. Sathyamurthy, J. Chem. Soc., Faraday Trans., 93, 773 (1997).
 (28) J.P. Gauyacq, J. Chem. Phys. 93, 384 (1990).
 (29) A. Bringer, J. Harris and J. W. Gadzuk, J. Phys.: Condens. Matter 5, 5141 (1993).
 (30) P. Hyldgaard, S. Hershfield, J. H. Davies and J. W. Wilkins, Annals of Physics 236, 1 (1994).
 (31) K. Flensberg, Phys. Rev. B 68, 205323 (2003).
 (32) A. Mitra, I. Aleiner and A. J. Millis, Phys. Rev. B 69, 245302 (2004).
 (33) Y. Imry, O. EntinWohlman, and A. Aharony, Europhys. Lett. 72, 263 (2005).
 (34) H. Ness, J. Phys.: Condens. Matter 18, 6307 (2006).
 (35) T. Markussen, R. Rurali, M. Brandbyge , and A.P. Jauho, Phys. Rev. B 74, 245313 (2006).
 (36) S. Roche, J. Jiang, F. Triozon, R. Saito, Phys. Rev. Lett. 95, 076803 (2005).
 (37) J.C. Lanczos, J. Res. Natl. Bur. Stand. 45, 225 (1950).
 (38) J. K. Cullum, R.A. Willoughby, Lanczos algorithms for large eigenvalue computations, SIAM (2002).
 (39) A. G. Borisov, S. V. Shabanov, J. Comput. Phys. 209, 643 (2005).
 (40) J. Sjakste, A. G. Borisov, J.P. Gauyacq, A.K. Kazansky, J. Phys. B 37, 1593 (2004).
 (41) D. J. Tannor, Introduction to quantum mechanics: a time dependent perspective, University Science, Sausalito (2006).
 (42) T. Holstein, Ann. Physs. (NY) 8, 325 (1959).
 (43) Y. Meir, N. S. Wingreen, Phys. Rev. Lett. 68, 2512 (1992).
 (44) P. K. Hansma, Phys. Rep. 30, 145 (1977).
 (45) B. Y. Gelfand, S. SchmittRink, and A. F. J. Levi, Phys. Rev. Lett. 62, 1683 (1989).
 (46) J. A. Stovneng, E. H. Hauge, P. Lipavsky and V. Spicka, Phys. Rev. B 44, 13595 (1991).
 (47) N. S. Wingreen, K. W. Jacobsen and J. W. Wilkins, Phys. Rev. Lett. 61, 1396 (1988); Phys. Rev. B 40, 11834 (1989).
 (48) J. W. Gadzuk, Phys. Rev. B 44, 13466 (1991).
 (49) J. Bonca, S. A. Trugman, Phys. Rev. Lett. 75, 2566 (1995).
 (50) E. G. Emberly and G. Kirczenow, Phys. Rev. B 61, 5740 (2000).
 (51) M. Cizek, M. Thoss and W. Domcke, Phys. Rev. B 70, 125406 (2004).
 (52) Ph. Durand, I. Paidarová and F.X. Gadea, J. Phys. B: At. Mol. Opt. Phys. 34, 1953 (2001).
 (53) A. Baratoff and B. N. J. Persson, J. Vac. Sci. Tech. 6, 331 (1988).
 (54) J.P. Gauyacq, Dynamics of negative ions, World Scientific Lecture Notes in Physics, Singapore (1987).
 (55) M. Persson and B. Hellsing, Phys. Rev. Lett. 49, 662 (1982); B. Hellsing and M. Persson, Phys. Scr. 29, 360 (1984).
 (56) We assume that the main difference between the negative PES and the neutral one is the change in geometrical conformation of the negative ion with respect to the neutral molecule. However, frequencies also change and as a matter of fact the electronvibration coupling should be computed for the negative PES.
 (57) J. Repp, G. Meyer, S. Paavilainen, F. E. Olsson and M. Persson, Phys. Rev. Lett. 95, 225503 (2005).
 (58) G. D. Mahan, Manyparticle physics, Kluwer Academic/Plenum Publishers, 3 Ed., New York (2000), p. 226.
 (59) T. Frederiksen, N. Lorente, M. Paulsson, M. Brandbyge, Phys. Rev. B 75, 235441 (2007).
 (60) N. A. Pradhan, N. Liu and W. Ho, J. Phys. Chem. B 109, 8513 (2005).
 (61) N. Ogawa, G. Mikaelian and W. Ho, Phys. Rev. Lett. 98, 166103 (2007).
 (62) J. R. Hahn, H. J. Lee and W. Ho, Phys. Rev. Lett. 85, 1914 (2000).
 (63) L. J. Lauhon and W. Ho, Phys. Rev. B 60, R8525 (1999).
 (64) N. Lorente, R. Rurali and H. Tang, J. Phys.: Condens. Matt. 17, S1049 (2005).
 (65) N. Lorente and M. Persson, Faraday Discussions 117, 277 (2000).
 (66) F. E. Olsson, N. Lorente and M. Persson, Surf. Sci. 522, L27 (2003).
 (67) L. C. Davis, Phys. Rev. B 2, 1714 (1970).
 (68) B. N. J. Persson and A. Baratoff, Phys. Rev. Lett. 59, 339 (1987).
 (69) M. Paulsson, T. Frederiksen, and M. Brandbyge, Phys. Rev. B 72, 201101(R) (2005).
 (70) O. Tal, M. Krieger, B. Leering, and J. M. van Ruitenbeek, condmatt.meshal arXiv:0801.3031v1
 (71) T. Frederiksen, M. Brandbyge, A. P. Jauho, N. Lorente, J. Comp. Elec. 3, 423 (2004).
 (72) C. Zhang, E. Jeckelmann and S. R. White, Phys. Rev. Lett. 80, 2661 (1998); Phys. Rev. B 69, 14092 (1999).
 (73) M. Nest and T. Klamroth, Phys. Rev. A 72, 012710 (2005).
 (74) P. Krause, T. Klamroth, and P. Saalfrank, J. Chem. Phys. 127, 034107 (2007).