Quantum dissipation driven by electron transfer within a single molecule investigated with atomic force microscopy

Intramolecular charge transfer processes play an important role in many biological, chemical and physical processes including photosynthesis, redox chemical reactions and electron transfer in molecular electronics. These charge transfer processes are frequently influenced by the dynamics of their molecular or atomic environments, and they are accompanied with energy dissipation into this environment. The detailed understanding of such processes is fundamental for their control and possible exploitation in future technological applications. Most of the experimental studies of the intramolecular charge transfer processes so far have been carried out using time-resolved optical spectroscopies on large molecular ensembles. This hampers detailed understanding of the charge transfer on the single molecular level. Here we build upon the recent progress in scanning probe microscopy, and demonstrate the control of mixed valence state. We report observation of single electron transfer between two ferrocene redox centers within a single molecule and the detection of energy dissipation associated with the single electron transfer.


Supplementary Note 1: Additional information about the formation of charge states and their characteristics
In this section, we provide additional information about the formation/annihilation of the single electron charging processes depending on an applied voltage and a tip position with respect to the molecule and limited possibility to discriminate the charge states in far tip-sample distances.
First, we describe an experiment, which demonstrates how the +1h charge state can be annihilated by a proper choice of tip position and applied bias voltage. Supplementary Figure 2 shows a sequence of constant height frequency shift and dissipation images acquired for different voltages. At bias voltage −2 V the molecule is in the +1h charge state, which is revealed by the presence of the characteristic line, see Supplementary Fig. 2a,d. In the next image, the bias voltage was reduced to a less negative value of −1.5 V while the tip height was kept constant. When the probe approaches the molecule, at one moment, an electron tunnels from the tip to the molecule and the hole is annihilated as shown on Supplementary Fig. 2b,e. Consequently, the sharp line disappears and the frequency shift becomes more negative, revealing large attractive force acting between tip and sample.
For the measurement shown on Supplementary Fig. 3, we first acquired a parabola in far distance (black line), where electron tunneling from tip to the molecule is not possible. Consequently, the charging process did not take place. After that, we approached the tip towards to the molecule by 1 nm and we swept the bias voltage from +4 to −4 V. In this case, we clearly observed two characteristic jumps in the Kelvin parabola revealing two charging events (see blue line). Next, we kept the bias voltage set to −4 V to maintain the molecule at +2h charge state and we retracted the tip to the initial tip sample distance by 1 nm. Finally, we acquired again the Kelvin parabola in this tip-sample distance (dashed red curve). We can see that the original and final Kelvin parabolas acquired in far distance are very similar and we cannot distinguish the presence of charge molecular states in such distances. We attributed this effect to the presence of strongly polar NaCl substrate, which may efficiently screen the local point charges placed on the surface. Unfortunately, signal-to-noise ratio of our instrument does not allow to detect variation of the frequency shift in orders of 0.01 Hz, which would reveal the +1h charge state in such far distances. This makes distinguishing the charge state of the molecule very difficult in far distances.
The lack of resolution in the far distance is actually already obvious from the vertical cut of frequency shift displayed on Fig. 2g of the main text. It is evident that the frequency shift signal fades away quickly when increasing the tip-sample distance. Thus characteristic round features indicating the position of the two ferrocene units are only visible in close distances, as one can see from Fig. 2g (x-y constant height and vertical x-z ∆f images, where you can see bisFc molecule and singularity line). In tip-sample distances above the dissipation line, the frequency shift signal is below the noise-level and we cannot make any statement about the charge state of the molecule.
In far distance, we are no more sensitive to the presence of charge states, as discussed previously. The variation of the frequency shift contrast at different tip heights can be seen from Supplementary Fig. 4. First we scan the molecule in neutral state (U = −0.1 V) in close distance, after that we switch the molecule to +1h charge state (U = −2 V) but keeping the tip height constant. The presence of the +1h state is clearly visible by the presence of the sharp line in both channels. After that, we retracted the tip 0.8 nm away and we lowered the bias voltage to 0.1 V, which in close distances would quench immediately the molecule to the neutral state. Nevertheless, the large tip-sample distance prevents the electron tunneling from the tip to the molecule and therefore the molecule remains in the +1h charge state. This is confirmed by the presence of the characteristic sharp line. However, we no more see the circular features representing ferrocene units in the frequency shift and the signal becomes monotonous except for the sharp line. Finally, we decreased the tip-sample distance and the molecule becomes neutral again. The fact that we can also modulate the position of the charge within the molecule-as revealed by the presence of the sharp line-even at low bias demonstrates that the electrostatic interaction is mainly driven by the presence of a static charge on the tip, but not by an induced one.

Supplementary Note 2: Description of the numerical model
The oscillating tip of the atomic force microscope indirectly measures forces between the tip and the sample through the shift ∆f of resonant frequency fr at which the cantilever with the tip oscillates. In addition, energy dissipated in the course of the oscillation can be detected and used to characterize (the non-conservative component of) the forces. Assuming a stiff cantilever, so that the tip motion, represented by time varying tip height zt, is only very weakly perturbed by the external forces and can be well approximated by the frequency shift ∆ is related to the force as ∆ = − r 2 ∫ ( )cos(2π r ) 1 r 0 (2) and the dissipated energy diss is diss = 2π r ∫ ( ) sin(2π r ) . 1 r 0 (3) In the above, A denotes the oscillation amplitude, fr the resonant frequency and F(t) the timedependent force. It is only a matter of convention that we chose the cosine function to describe the tip oscillation, rather than e.g. the sine function or a sine function with some general phase-shift. This choice just fixes the origin of the time scale, t = 0. In the following, we will only be interested in the components of F(t) which are related to the extra charge localized on the molecule. All other components of the force are typically conservative, hence not contributing to Ediss, and only contribute a smooth featureless background to ∆f. The charge-dependent component of the force where FA (FB) denote the electrostatic force felt by the tip because of the extra charge +1h at location A (B) in the molecule; pA and pB are the time-dependent probabilities that the respective locations A, B are occupied by the charge, and the coordinates xt, yt, zt specify the position of the tip.
If we represent the electrostatic effects of the tip by a single point charge Qt placed at (xt,yt,zt), the forces FA and FB are given by Note that only the z component of the Coulombic force is relevant for the motion of the tip. The electrostatic interaction of such tip then shifts the energies that correspond to the charge localized at sites A or B. These energy shifts are again dictated by the Coulomb interaction, so respectively.
When the charge on the molecule switches from one location to another, the molecule reacts by reorganizing its structure. The reorganization in general involves both a shift in positions of atomic nuclei as well as a redistribution of electron clouds which takes place in order to partially screen the localized charge. Let us assume that the slowest component of the reorganization occurs at a time scale comparable to the tip oscillation and that it can be characterized by a single reaction coordinate, which we are going to denote q; see Fig. 3a,b in the main text. The energy landscape in q will be assumed to be described by a parabolic potential VA(q) or VB(q), the minimum of which depends on whether the molecular charge resides at site A or B: where δ is the distance in the q between the equilibria that correspond to sites A and B, respectively, being occupied; mq is the effective mass and ωq the angular frequency of the oscillator formed by the potential wells VAB(q). We used the freedom to arbitrarily define the origin of the q axis in order to place the potentials VA and VB in a symmetric way with respect to the origin. We are going to characterize the statistical state of the molecule by probability densities ρA(q,t) and ρB(q,t), The q-dependent energy change accompanying the switch of the charge location from B to A is where ∆E0 is the initial energy difference between the two possible charge positions A and B (with the molecular structure in the corresponding relaxed configuration for both of these charge positions) in the absence of the tip. A non-zero initial energy difference ∆E0 is a consequence of a slight asymmetry of the molecule adsorbed on the insulating surface. Such an asymmetry is necessarily present because the distance of the Fe redox centers in the molecule is incommensurate with respect to the NaCl(100) surface structure, so the two parts of the molecule with the respective redox centers will adsorb in a slightly different geometry and they will experience different electrostatic fields and dielectric screening from the ionic substrate. The immediate energy difference ∆EAB(q,t) is time dependent through its dependence on tip position, in particular the vertical tip coordinate zt(t). Note that the quantity Δ ( ), which appears in equation (2) of the main paper, collects those terms of ∆EAB(q,t) that depend on the position of the redox centers with respect to both the tip and the substrate: Δ ( ) = A − B + Δ 0 .
We model the dynamics of the probability densities in the potential wells by the Smoluchowski diffusion equation in order to mimic the relaxation of the molecule towards its equilibrium configuration for a given charge location. Additionally, q-preserving charge hopping from A to B and back will be allowed with rates kBA and kAB [as in Valkunas, L., Chmeliov, J., Krüger, T. P. J., Ilioaia, C. & Grondelle, R. JPC Lett. 3, 2779-2784 (2012)]. The resulting coupled differential equations for ρA(q,t) and ρB(q,t) (adapted from Mukamel: Principles of Nonlinear Optical Spectroscopy book, Eq. 12.48) become where kB is the Boltzmann constant, T is the thermodynamic temperature of the substrate with which the molecule equilibrates, γ is the friction constant for the reorganization coordinate, and the hopping rates are Supplementary Equation (12b) corresponds to equation (2) from the main paper, while (12a) gives the hopping rate in the opposite direction. The parameter Λ in the above equations is the part of reorganization energy that pertains to degrees of freedom (DOF) which relax much faster than those represented by coordinate q. These fast DOF are those which we dubbed "Marcus" DOF, while those represented by q are the "explicit" DOF. The coefficient k0 scales with temperature and Λ according to The notation in Supplementary Eqs. (11) can be simplified if we introduce a characteristic time constant KΛ for the diffusion and a temperature-dependent diffusion parameter D(T) by and The coupled dynamical equations for the probability densities then become The number of independent parameters in the above equations can be further reduced by rescaling q according to q̃ = q/δ and introducing a new parameter The probability densities shall be rescaled accordingly, ρÃ(q,t ) = ρA(q,t)δ and ρB(q,t ) = ρB(q,t)δ.

Supplementary Eqs. (16) can now be rewritten as
and the expressions for the potential wells, Supplementary Eqs. (8), as Although we have omitted the tildes in the above Supplementary Eqs. (18) and (19), it was already written in the rescaled quantities. We may interpret the new parameter EΛ as a reorganization energy associated with the relatively slow dynamics in coordinate q (in contrast to the reorganization energy Λ, associated with the "faster" reaction coordinates). The dynamics of the electron transfer (ET) is then essentially characterized by four parameters (apart from the temperature, which we assume to be known from measurement): the diffusion time scale KΛ, the slow reorganization energy EΛ, the fast reorganization energy Λ, and the hopping-rate coefficient k0 (or, equivalently, the transfer matrix element M).

Supplementary Note 3: Sensitivity of the line shape of ∆f(z), Ediss(z) signals to the electrostatic model of the tip
The shape and spatial location of the sharp linear feature related to the electron transfer is dictated by the interaction of the charged centers in the molecule with the tip of the microscopic sensor and with the substrate. Electrostatic properties of the tip termination are the decisive factor here. Supplementary Figure 5 demonstrates the dependence of the line shape on the tip by showing how the shape of the imaged feature changes after the tip is intentionally modified by bringing it to a close contact with the substrate and the same molecule is then rescanned with the modified tip. Supplementary Figure 6 shows the same effect for the simulated data. Asymmetric tips are modelled with two-point charges on the apex instead of just one.

Supplementary Note 4: Numerical details regarding simulation of ET dynamics
According to Supplementary Eqs. (12), ET is most probable to occur when ΔEAB(q, t) = ± (the sign depending on the direction of the transfer). In such a situation, the energy dissipated instantaneously (into the "Marcus" DOF) as the charge hops from one site to another equals the reorganization energy Λ. Another contribution to the dissipation of energy Ediss is the damped relaxation to the ground state in the parabolic potentials (the explicit DOF). This contribution scales up with the other reorganization energy EΛ and also increases with the diffusion time constant KΛ as the relaxation gets faster with larger KΛ.
We have carefully chosen the parameters of the numerical model in order to reproduce the experimental conditions and the observed signal as faithfully as possible. The resonance frequency and stiffness of the sensor, fr = 1 MHz, k = 1 MNm -1 , were dictated by the known properties of the measuring device. All simulations presented here were done for an amplitude of A = 40 pm, for which we achieved the best sensitivity in the experiment (an optimal trade-off between a signal to noise ratio, which improves with increasing amplitude, and a loss of spatial sensitivity for larger amplitudes). The electric charge on the tip apex, Qt = +0.25 e (where e is the elementary charge) was chosen in order to roughly reproduce the typical magnitude of the measured signal on the sharp features related to charge hopping (frequency shift in the order of 0.1 Hz, dissipated energy several meV). The initial (tip-independent) site asymmetry of ΔE0 = 10 −20 J ≈ 62.4 meV then nicely reproduces the curved shape of the linear feature observed in some of the experiments; see Given the choices made so far, the parameters which directly relate to the dynamics of the electron transfer were adjusted with the following conditions in mind: (i) The effective ET rate <kET> derived from the experimentally-determined ratio of the frequency shift signal to the dissipation signal according to equation (1) of the main paper fell, as Fig. 2k testifies, into the range from 2×10 6 s -1 up to about 10 7 s -1 , with the exception of one data set that went up to about 1.8×10 7 s -1 . We tried to get the same order of magnitude for the corresponding quantity derived from the simulated signal. (ii) The results should be rather insensitive to temperature; refer again to Fig. 2k to see this is the case for the measured values when changing the temperature from about 5 K to about 1 K. (iii) The shape of peak profiles in z-scans as seen in both signal channels, frequency shift as well as dissipation, should be a bell-like curve resembling a Gaussian or Lorentzian function (these two are not distinguishable within experimental precision) rather than being almost rectangular or semi-circle shaped, as it could happen with an inappropriate parameter choice in our model. Although we were not able to make a uniquely optimal choice of the parameters, here we give the values that fulfilled the above criteria reasonably well. The reorganization energy for the explicit degrees of freedom EΛ = 4×10 −23 J ≈ 2.5 meV, relaxation rate of the explicit DOF KΛ = 2×10 6 s −1 , Marcus reorganization energy Λ = 4×10 −22 J ≈ 2.5 meV, the transfer rate k0 = 7.07×10 7 s −1 for the temperature of T = 1.2 K and k0 = 3.46×10 7 s −1 for the temperature of T = 5 K, which corresponds to the transfer matrix element M = 3.65 μeV. These values are adopted to obtain the results shown in Figures 3 and 4 in the main text.
The above choice of parameters KΛ and EΛ for the ET dynamics corresponds to a reorganization intermediate rate of the explicit DOF relaxation. In this case, the relaxation of the explicit DOF takes place on a time scale comparable to the oscillation period of the AFM tip. Let us now compare our full model with the above parameters (at 5 K) to a simplified version that is based solely on the Marcus model for charge hopping while it completely disregards the explicit degree of freedom. Within such model, the probabilities that the charge is located at a particular position evolve according to  (20), but the shape of the ∆f, Ediss peaks is more rounded. In Supplementary Figure 11a,b, we compare the results from the full model to the simplified all-Marcus model. The reorganization energy for the all-Marcus model was chosen to be Λ = 6×10 -22 J ≈ 3.75 meV and the transfer rate was set to k0 = 2.45×10 7 s -1 (at T = 5 K). Note that the value of Λ was chosen to lie in between Λ and Λ + EΛ of the full model. Such a choice ensures an approximate agreement of the dissipated energy Ediss between the two models: The typical dissipation energy Ediss as a result of one charge-transfer event in the Marcus model is roughly the reorganization energy Λ. In the full model, on the other hand, the typical dissipated energy Ediss equals the complete reorganization energy of the fast (Marcus) DOF (Λ) plus part of the reorganization energy of the slower explicit DOF (EΛ). The reason why only part of the reorganization energy of the explicit DOF dissipates is that these slower DOF do not have the time to fully reorganize before next electron transfer occurs. In spite of a crude agreement of signal magnitudes, peak shapes Ediss clearly differ between the full and the simplified Marcus model. The full model with the explicit DOF reproduces the shape of the Ediss peaks in experimentally detected signal somewhat better.
Yet more prominent discrepancies in peak shapes arise with respect to the experiment if we used another, extremely simplified ET model. This model consists in the assumption that charge the immediately switches from one site A to B as soon as a difference in energies corresponding to these sites, Δ ( ) = A − B + Δ 0 , surpasses a certain constant threshold δE. In the same way, the charge moves from B to A as soon as Δ ( ) < − . When −δ < Δ ( ) < +δ , no electron transfer can happen. Simulated ∆f(z), Ediss(z) profiles and their comparison to the full model are shown in Supplementary Fig. 11c,d; choosing δE = 4×10 −22 J ≈ 2.5 meV as the threshold energy. In this threshold model, the dissipated energy Ediss in one cycle is always either 2δE or zero. Thus, the dissipation signal Ediss appears only when the oscillating tip modulates the energy difference ∆E(t) between the two redox sites larger than the threshold energy δE during tip approach and retraction and consequently the electron transfer can proceed in both directions. On the contrary, the zero-dissipation case occurs if the tip-induced energy variation ∆E(t) does not suffice to drive the electron transfer. This on/off character results in rectangular shapes of the peaks in the dissipation Ediss channel of the measured signal, as shown on Supplementary Fig. 11d. So on the one hand, by comparing simulated ∆f(z) profiles for three different ET models, see Supplementary Figs. 11a,c, we can realize that the shape of ∆f(z) is rather insensitive to the ET models. On the other hand, it is evident that the shape of Ediss peaks is sensitive to the ET dynamics. Thus, it can provide more insight into the ET dynamics. In order to relate the shape of the Ediss(z) peaks to ET rates, we plot in Supplementary Fig. 12 the mean number of ET events that take place during one oscillation cycle of AFM probe for each model. Variability of transfer probability peak shapes in these plots correlates very well with the variations of the dissipated energy Ediss(z) profiles. This demonstrates the direct relation between the ET dynamics and the shape of Ediss(z) profiles.

Supplementary Note 5: Additional discussion on the explicit and Marcus DOF
A better agreement of our full numerical model with the experimental evidence as compare to what we call the Marcus model, reveals the importance of the explicit DOF to describe the electron transfer process. The better agreement concerns the calculated line shapes and their amplitude dependence. On the other hand, the very existence of energy dissipation and its weak temperature dependence could have been explained by the Marcus model as well.
In principle, there are two components which define the out-of-equilibrium state of our system. They can be identified in the time evolution of the wave packets shown in Figure 3 in the manuscript.
The first component is related to a non-zero occupation probability pA or pB of a site that is higher in energy. This component is also included in the Marcus model and it is related to dynamical motion of the probe and to a finite value of the transfer rates kAB and kBA. Here, in principle, the time scale of the DOF, into which the energy dissipates simultaneously with the charge transfer, is faster than the probe oscillation. The dissipation mechanism is related to irreversible energy transfer from the molecule to the substrate. We expect that this fast part comprises mostly of the intramolecular vibrational modes of the studied molecule, the near equilibrium motion (phonons) of the substrate, etc.
The second component, which defines the out-of-equilibrium condition, consists of a location of the center of mass of the wave packets out of the minimum of the parabolas. This component is only present in the full model. Importantly, this explicit DOF additionally modifies the occupation probability. Consequently, it also influences the line shape of the frequency shift and the dissipation signal, which allows to improve the agreement with the experiment. The comparison of the numerical model to experiment suggests the presence of some slow component of the bath dynamics. We do not have any definite physical explanation regarding its origin. It could be related, for instance, to the change of adsorption site of the molecule on NaCl or to a complex structural reorganization of the ionic substrate. Actually, the latter possibility is supported by the DFT calculations of charged molecules, showing pronounced relaxations of the substrate (bond length variation larger than 0.1 Å). In general, there may be other degrees of freedom than just the electron location that are driven out of equilibrium by the tip-oscillation driven charge transfer.

Supplementary Note 6: Other possible numerical models of the system
Finally, we would like to note that the numerical model we presented in this work may not be the unique way how to describe the weakly dependent ET process in our system. For example, the transfer rate as derived by J. Jortner [J. Chem. Phys. 64, 4860-4867 (1976)] is temperature independent in the low temperature limit. However, in this model all modes are effectively frozen and there is no dissipation. We have considered Jortner theory in an early stage of our theoretical analysis of the experiment, but the lack of dissipation excludes it from being the right answer to our particular problem. In our particular case, the explanation does not require temperature independent rates, rather the discrepancy with respect to equilibrium predictions arises from nonequilibrium (driven) state of our system. It is important to note, that if we assume that all the components of the environment relax fast on the time scale of the electron transfer process, i.e. we can apply Marcus theory to all degrees of freedom of the environment, the calculation still gives reasonable results. We can see however, that a more general theory, which allows for some degrees of freedom of the environment to relax relatively slowly, gives better agreement with the measured and amplitude dependence.
In principle, other formalisms designed to describe the dynamics of the open quantum systems such as Caldeira-Leggett could be applied to the description of the experiment. We hope the presented experimental data may stimulate further research activities applying methods describing quantum open systems to this system.

Supplementary Note 7: Background subtraction of ∆f(z) in the numerical model
In order to compare experimental and simulated data, it is essential to subtract the background on which the frequency dips are superimposed in the ∆f(z) spectra, because the sources of the background differ between the experiment and the simulations. We have described how we processed the experimental spectra in the Methods section of our paper. Here, we discuss the issue of the background subtraction for the simulated ∆f(z) spectra. Our numerical model does not contain any long-range van der Waals forces, so we do not need to subtract them. However, there is still a force originating from the point charge sitting on the molecule. Only the periodically tipdriven hopping of this charge between two redox centers causes the peaks, while the frequency shift related to the z-gradient of a static force that originates from the charge sitting on one particular site should be considered part of the background.
In the next, we describe a procedure to get rid of this background component of the simulated frequency shift. From the simulation, we determined time averages of occupation probabilities pA(t) and pB(t). This was done separately for each point of the spectrum, that means for each value of height z0 around which the tip oscillated. Afterwards, we modelled the frequency shift again, imposing a constant fractional charges on site A and B, which correspond to the time-averaged values of pA and pB, respectively. We took this frequency shift ∆f(z) profile to be our background and we subtracted it from the frequency shift ∆f(z) given by our original simulation. The reader can convince oneself in Fig. 3c,d or Supplementary Fig. 9 that the simulated spectra consist exclusively of the peaks surrounded by pure zero signal, thus comparing well to the backgroundsubtracted curves from the experiment.

Supplementary Note 8: Relation between the effective transfer rate <kET> and time dependent rates kAB, kBA
Our definition of the effective transfer rate <kET>, Eq. (1), follows analogous expressions that were independently suggested by the group of P. Grutter [Roy-Gobeil, A., Miyahara, Y. & Grutter, P. Nano Lett. 15, 2324-2328(2015], their Eq. (6), and ourselves [Ondráček, M., Hapala, P. & Jelínek, P. Nanotechnology 27, 274005 (2016)]. The applicability of this equation is limited to cases with a pronounced dissipation. This condition is only fulfilled when the switching rate is similar to the oscillation frequency fo of the probe. In other words, the dissipation signal becomes maximized when switching occurs only a few times during one oscillation cycle and when the time-dependence of the occupancy of two possible states is in phase delay with respect to the oscillation amplitude.
The effective rate <kET> relates to the phase shift between the tip oscillation and the electrostatic force FQ(t) sensed by the probe. Tip oscillation directly translates into periodic variations of the energy level difference Δ ( ) = ( , , ) − ( , , ), while FQ(t) is given by the immediate location of the hopping charge. Continuous variation of the energy level difference ΔE(t) results into temporal variation of the charge occupancies on sites A and B, i.e. pA(t), pB(t), but with a certain time delay with respect to the oscillation amplitude of the probe due to finite charge transfer rates kAB and kBA. This delay causes a phase shift between the oscillating tip height zt(t) and the force FQ(t), as the latter depends on the charge occupation probabilities pA(t) and pB(t).
The link of the effective switching rate <kET> to the phase shift can be demonstrated by looking at the Fourier component of FQ(t) that corresponds to the base cantilever frequency fr. Suppose we take Q ( ) ∝ cos (2 r + ), substitute this into the expressions for the frequency shift and the dissipated energy, equations (2) and (3) of this SOM, respectively. We interpret the resulting Δf and Ediss as the peak values associated with the charge switching signal and insert them into the definition of <kET>, Eq. (1) from the main paper. We obtain < > = r cotg ≈ 1 2 ⁄ ; where φ is the phase shift with respect to the tip position zt(t) and therefore also with respect to ΔE(t), while τ is the mean time it takes for the charge to switch into the new location after the tip passes the position where ΔE = 0. The condition ΔE = 0 defines the situation in which the bottom of both parabolic potentials for the two redox sites are aligned on the same level, as depicted in Fig. 3b of the main text. The approximate expression for <kET> that involves the mean delay τ is valid as long as the charge position switches in (almost) every cycle of the tip oscillation and the switching happens relatively soon after ΔE = 0 as compared to the whole oscillation period. Such conditions correspond to a small phase shift φ. That the assumption of small φ was indeed fulfilled could be checked a posteriori by observing that <kET> evaluated from our measured data was several times larger than fr. Also note that our <kET> corresponds to one half of ΓΣ from [Roy-Gobeil, A., Miyahara, Y. & Grutter, P. Nano Lett. 15, 2324-2328(2015], eq. (6). Finally, let us note that <kET> would have a straightforward interpretation in terms of the immediate transfer rates kAB and kBA only if the sum kAB + kBA were constant (or at least only slowly varying with tip oscillations), in which case <kET> ≈ (kAB + kBA) / 2. This however definitely does not hold in our model of the transfer mechanism, as we discuss further below.
The rate <kET> is a characteristic time scale of the charge switching, one which can be extracted from the experimental data in a straightforward manner. However, the interpretation of <kET> in terms of the physical parameters of our theoretical model is far from being straightforward. This in particular holds about the connection between <kET> and the actual time-and configurationdependent switching rates kAB(q, t) and kBA(q, t), as given by Eq.
(2) in the main paper or Supplementary Eqs. (12). On the one hand, the rates kAB and kBA go from 0 to ~20 MHz in our model, as one can see in Fig. 3c, while the effective rate <kET> calculated from Eq. (1) gives values in the range of 6-8 MHz for different tip-sample distances, so these quantities are at least of the same order of magnitude. On the other hand, the dependence of <kET> on kAB(q, t) and kBA(q, t) is significantly nonlinear.
The crucial observation for establishing a relation between these rate quantities is the inverse proportionality between <kET> and the mean delay τ of the charge switching. So as not to complicate the discussion any further, we will assume here that, in accord with Supplementary Eqs. (12) in Supplementary Note 2, the rate kBA is the same function of ΔEAB as the opposite rate kAB is of minus ΔEAB. Furthermore, most of the time, switching in one direction is much more probable than in the opposite direction. So let us only care about one of the possible directions of the switching when estimating τ, say from B to A. The DOF described by the coordinate q is not directly related to the force FQ(t). Only the location of the charge either on site A or site B decides about the force felt by the tip. We may therefore consider rates integrated over q, such as kAB(t) = ∫dq ρB(q) kAB(q, t), for the present discussion. Note that this is also the meaning of the rates shown in the middle panel of Fig. 3c in the main paper.
Now, there are two interwoven trends how kAB(t) affects τ and thus <kET>. First, <kET> in general grows with growing kAB(t): when kAB is larger, the transfer of the charge to the other redox center occurs sooner, so τ becomes shorter and <kET> larger. If kAB(t) were roughly constant during the relevant phases of tip oscillation in which switching from B to A occurs, <kET> ≈ kAB(t) would hold. Second, in our model, kAB(t) is in fact not constant but rather changes, often quite abruptly, and peaks for some ΔE = ΔEpeak, somewhere in the range < |ΔEpeak| < ( + E ). If the peak of the transfer rate is sharp (typically when kBT << ), we may assume that the charge switches exactly when ΔE = ΔEpeak and then <kET> ≈ 2πfr (∂E/∂z) A/ΔEpeak, where A is the oscillation amplitude. If, however, the rates kAB (kBA) exhibit rather broad peaks as functions of ΔE, as it is also the case in our study, the two above described trends participate in determining the value of the effective rate <kET> in a rather tangled way.
To conclude, in our model the effective rate <kET> expresses, in a certain way, the phase shift between tip oscillation and the electron transfer switching. Although the phase shift is associated with the kAB(t) and kBA(t) rates, <kET> depends not only on the magnitude of kAB(t) and kBA(t), but also on how they change in time during tip oscillation. The sooner kAB(t) or kBA(t) grows after passing the equilibrium position, the larger <kET> becomes. Moreover, the effective rate <kET> depends non-linearly on how large kAB(t) and kBA(t) rates are. In other words, <kET> gives only an approximate estimate of the switching rate in our model and it cannot provide direct information about time-dependent character of the kAB(t) and kBA(t) rates employed in the model.

Supplementary Note 9: Amplitude dependence of ∆f(z) and Ediss(z) in experiment and theory
To see how the character of singularity line depends on the oscillation amplitude, we performed series of constant-height measurements, where the oscillation amplitude was gradually increased. Supplementary Figures 13 and 14 show variation of the frequency shift and the dissipation signal in the amplitude range of 40-90 pm and 20-80 pm, respectively. In these amplitude ranges, we observe a modest dependence of the signal on the amplitude. Namely, the width of the frequency shift and the dissipation peaks increases slightly with the amplitude. In the case of the dissipation signal, an enhancement of the maximum value of dissipation signal with the amplitude is more pronounced than in the case of the frequency shift which remains almost constant. Note, tip height is not exactly maintained during variation of the oscillation amplitude due to a feedback response. Thus, the images shown on Supplementary Fig. 13 and 14 are not taken with exactly the same tip height.
In the next, let us compare the numerical results obtained with the full and Marcus model. Despite the fact that the theoretical plots show the tip-height profile, while the experimental profiles are displayed as lateral profiles acquired at the constant tip-height, we can directly compare the shape of the peaks (mainly width and the maximal value of the peaks) with respect to the oscillation amplitude.
In both models, we observe a gradual broadening of the peaks in both channels and increase of the dissipation peaks as the oscillation amplitude increases, in good agreement with the experimental evidence. However, the full model with both the Marcus and explicit DOF shows better match with the experimental evidence. Namely, the experimental data presented in Supplementary Fig.  13 and 14 are taken in far tip-sample distance regime, as the characteristic circular features reflecting the ferrocene units are not visible. The evolution of simulated frequency shift peaks with the amplitude in far tip-sample distances (two peaks on the right in Supplementary Fig. 11a) matches perfectly the corresponding evolution observed in experiment. First, the maximum of the frequency shift increases from A = 20 pm to 40 pm, while for the amplitude A = 80 pm it drops again and the peak becomes broader. This behavior cannot be reproduced in the Marcus model, where the frequency shift monotonically decreases with larger amplitude, see Fig. 5c. We also think that the evolution of the dissipation signal with the oscillation amplitude is better captured by the full model. Thus, the amplitude dependence provides another indication, that the full model with the Marcus and explicit DOF describes better the electron transfer mechanism, which takes place in our system.

Supplementary Figures
Supplementary Figure 1. Register with substrate. Constant height frequency shift image of bisFc molecule with atomic resolution of NaCl substrate. Centers of circles denote Cl atoms. Image was taken in two different tip -sample heights. First, the resolution on substrate was taken, then the tip was retracted by 0.6 nm and the bisFc molecule was imaged. Scanning area is 4.5 × 4.5 nm 2 , oscillation amplitude A = 40 pm, bias voltage U = 0 V. Series of Δf(U) parabolas taken in 2 different tip-sample heights. Black solid line represents a parabola taken above bisFc molecule in the neutral state at relative height z =+1 nm. Blue solid line represents a charging parabola taken at relative height z = 0 nm (1 nm closer to the bisFc molecule); bisFc becomes twice positively charged during its acquisition. Red dotted parabola is taken at relative height z =+1nm above twice positively charged bisFc molecule. Oscillation amplitude is A = 50 pm. , bias corresponds to a neutral molecule) and U = −2 V (corresponds to a once positively charged molecule, panels (b), (f)). The scanning area is 2 × 2 nm 2 , oscillation amplitude A = 50 pm for each image, temperature 5 K. Scanning direction in all images is from top to bottom.

Supplementary
Step 1 (a,e) presents a neutral molecule at a relative tip-sample distance z = 0 nm.
Step 2 (b,f) shows a once positively charged molecule (bias voltage U = −2 V) at same relative tip sample distance z = 0 nm. Between step 2 and step 3 the tip is retracted by z = +0.8 nm while keeping bias voltage U = −2 V.
Step 3 (c,g) represents charged molecule at "far" distance, relative tip sample distance z = +0.8 nm, scanned at bias voltage U = −100 mV that corresponds to bias were the molecule is normally found in the neutral state, but in this particular case discharging is prevented by the tip-sample distance.