Multidimensional Hydrogen Tunneling in Supported Molecular Switches: The Role of Surface Interactions

The nuclear tunneling crossover temperature ($T_c$) of hydrogen transfer reactions in supported molecular-switch architectures can lie close to room temperature. This calls for the inclusion of nuclear quantum effects (NQE) in the calculation of reaction rates even at high temperatures. However, standard computations of NQE relying on parametrized dimensionality-reduced models, quickly become inadequate in these environments. In this letter, we study the paradigmatic molecular switch based on porphycene molecules adsorbed on metallic surfaces with full-dimensional calculations that combine density-functional theory for the electrons with the semi-classical ring-polymer instanton approximation for the nuclei. We show that the double intramolecular hydrogen transfer (DHT) rate can be enhanced by orders of magnitude due to surface fluctuations in the deep tunneling regime. We also explain the origin of an Arrhenius temperature-dependence of the rate below $T_c$ and why this dependence differs at different surfaces. We propose a simple model to rationalize the temperature dependence of DHT rates spanning diverse fcc [110] surfaces.

Nuclear tunneling is an inherently quantummechanical process that can strongly impact the properties of matter in a wide variety of situations, ranging from biological enzymes to organic-based technologies [1][2][3][4]. In complex environments, it has been recognized that signatures of tunneling on rate processes are often not well captured by textbook theories [5,6] and in particular, in hydrogen transfer reactions, the small mass of hydrogen makes tunneling pronounced [7]. Still, a theoretical description of nuclear tunneling that goes beyond a simple one-dimensional approximation and considers anharmonic coupling between many degrees of freedom in larger-scale systems remains a challenge [8][9][10]. To build up a systematic understanding of multidimensional rate processes in the deep tunneling regime and the impact of the environment on hydrogen dynamics, a fully ab initio treatment of simple, yet non-trivial, reactions in well controlled conditions is desired.
Tetrapyrrole macrocycles, like porphyrin, naphtalocianine, and porphycene, show a remarkable diversity of functional properties [11][12][13][14][15], and were proposed, among others, as molecular switches [16][17][18]. The intramolecular hydrogen transfer reaction that occurs in the inner cage of these molecules, known as tautomerization, can be triggered remotely by different external stimuli [19][20][21][22][23][24]. Moreover, because the reaction takes place without a pronounced conformational change, these molecules can be incorporated in nanoscale devices. Within this area, controlling the dynamical properties of these molecules is central to advance rational design. Molecular diffusion and rotations have been more often addressed [25][26][27][28][29], while key aspects of the hydrogen transfer reaction mechanism and its temperature dependence remain poorly understood. This is especially due to the challenges in the description of multidimensional quantum dynamics that cannot be captured by perturbative treatments of anharmonic couplings [5,[30][31][32][33].
In this Letter, we study the effect of tunneling on the double intramolecular hydrogen transfer (DHT) of two representative systems, namely porphycene on Cu(110) and Ag(110) surfaces. In these reactions, the tunneling crossover temperature, which represents the temperature below which tunneling becomes greater than classical hopping over the barrier (T c = ω b /2πk B , where ω b is the imaginary frequency of the unstable mode at the transition state geometry [34]) lies close to room temperature. Specifically, we i) identify the multi-dimensional DHT mechanisms at different temperatures, ii) clarify their temperature dependence, and iii) elucidate the role of surface fluctuations in the deep tunneling regime. The results we obtain are able to explain puzzling experimental measurements [4,21] that showed an unexpected temperature dependence of the DHT. In addition, they show that instead of acting only as a passive observer of the reaction, in certain situations the surface takes a prominent role in the tunneling event.
We here employ a combination of density-functional theory (DFT) for the electronic degrees of freedom with the ring polymer instanton (RPI) approximation [35,36]  for the nuclear degrees of freedom. RPI can be viewed as an extension of Eyring [37] transition state theory (TST) which includes tunnelling and captures anharmonic contributions along the reaction pathway. It is a semi-classical method that uses discretized closed Feynman path (CFP) integrals to evaluate tunneling rates in the deep tunneling regime. This approximation finds the dominant stationary pathways in the CFP that connect reactants and products, the instanton pathways. Then, a series of steepest descent approximations to evaluate the flux-side correlation function leads to the instanton approximation for thermal reaction rates k inst [38]. The rate can be expressed as where A inst accounts for the harmonic fluctuations around the instanton pathway x inst , S is the Euclidean action, and β = 1/k B T with k B the Boltzmann constant and T the temperature. In the discretized CFP space, x inst is a first order saddle point. Albeit approximate, this method shows the best tradeoff in situations where the quantum exponential wall would prevent a full dimensional evaluation of the exact tunneling rate [6,39].
Within the Born-Oppenheimer (BO) approximation, evaluating Eq. 1 requires the BO energies and forces. We here employ the most accurate level of theory affordable for the system size and number of system replicas required. We employ DFT with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [40] including the Tkatchenko-Scheffler [41] dispersion correction modified to treat physisorption on surfaces [42]. We thereby ensure a good description of the anisotropic electron density redistribution upon porphycene adsorption [43,44]. At this level of theory, we do not expect reaction barriers to be quantitatively accurate [45], but expect to capture qualitative trends and the correct physics. In the supporting information (SI), we report selected geometries calculated with a range-separated hybrid functional for comparison. These calculations were enabled by the combination of the FHI-aims [46] all-electron code and the i-PI [47,48] universal force engine. Details and convergence tests are provided in the SI.
In porphycene, the inner-cage hydrogen atoms can adopt different configurations, giving rise to two stable tautomers, coined cis and trans. Each tautomer is two-fold degenerate and their structures on Cu(110) and Ag(110) are shown in Fig. 1. On these surfaces, the cis tautomer is more stable than the trans tautomer by 172 and 20 meV, respectively. Indeed, on the same surfaces, experiments observe almost exclusively the cis tautomer [4,23,49]. The DHT of the cis→cis reaction can happen through two possible mechanisms. The concerted mechanism, where the hydrogen atoms are transferred together, without the existence of a stable intermediate and passing through a second-order saddle point, and the stepwise mechanism, where the hydrogen atoms are transferred sequentially and the reaction involves a trans intermediate.
We have recently studied the DHT of porphycene in the gas phase [33], and found that there is a competition between a concerted and a stepwise DHT at different temperatures. However, new effects arise due to the interaction between the molecule and the surface. In order to rationalize them, we divide the influence of the surface on the physisorbed molecule into static and dynamical effects. The static surface effects refer to the change in the potential energy landscape upon adsorption, and can affect both static and dynamical properties of the adsorbed molecule. For example, while the global minimum in the gas phase is the trans tautomer in a flat conformation, on Cu(110) and Ag(110) the relative tautomer stability is reversed and the molecule is buckled, as shown in Fig. 1b and 1c. This distortion changes the DHT energy barriers and consequently the hydrogen dynamics. Dynamical surface effects, on the other hand, refer to the impact that the motion of the surface atoms have on the molecular properties. As it will be shown, they can significantly influence the dynamics of the DHT.
On Cu(110), we calculate T c = 264 K for the DHT and show the instanton tunneling rates of the cis→cis reaction in Fig. 2a. We show rates between 75 and 150 K because they lie well within the deep tunneling regime and we wish to compare with experimental results in the same temperature range [21], shown in the inset. In this At these temperatures, kinst only includes contributions from the stepwise mechanism . The inset shows experimental results from Ref. [21]. (b) Calculated DHT kinst for the cis → cis reaction of porphycene on Ag(110) between 7.5 and 110 K. Orange circles represent kinst of the stepwise mechanism and blue squares of the concerted mechanism. The inset shows experimental results from Ref. [4]. (c) Schematic 1D potential energy surfaces of the stepwise (left) and concerted (right) reactions. The ZPE for reactant and products are shown, while the ZPE at the barrier top is not (but it was included in the multidimensional calculations). (d) Schematic representation of the different temperature dependence regimes of the stepwise (orange) and concerted (blue) DHT rate of adsorbed porphycene. Full/dashed lines represent the dominant/minor mechanism. The numbers 1, 2 and 3 illustrate the three different regimes, see text.
case, we find a very good agreement with experiment throughout the six-orders of magnitude variation of the rate. The calculated effective activation energy E A is 190 meV. This compares well with the value of 168 ± 12 meV reported by experiments. Both the calculated rates and the measured ones show an Arrhenius-like temperature dependence even though tunneling is dominant at this temperature range. We proceed to explain the origin of this dependence.
The calculated DHT rates contain contributions almost exclusively from the stepwise mechanism, because at these temperatures it is several orders of magnitude faster than the concerted alternative. We propose that the trans intermediate was not observed because its predicted residence time is ≈ 0.1 ns, which lies beyond the time resolution of STM experiments (≈ 100 µs). The instanton trajectory provides an intuitive view of the reactive process by showing the main "instantaneous" tunneling configuration that the delocalized nuclei adopt. The trajectory is visualized in Fig. S1 at several temperatures. As the temperature is decreased, the reaction takes shorter pathways and crosses regions of higher PES energy. Moreover, we observe a considerable contribution from heavier atoms like C and N and, interestingly, even Cu atoms to the tunneling mechanism (see Tab. S11). These are manifestations of the multidimensional nature of the tunnelling process and prove that reducing the problem dimensionality is inappropriate [8,50,51]. After these considerations, the observed E A can be understood as follows. At low enough temperatures, when only the vibrational ground states (VGS) are populated and a further decrease of the temperature does not affect the vibrational populations, the trans→cis reaction will proceed from the VGS of the reactant and will be constant with temperature. As a consequence and because of detailed balance, the inverse reaction cis→trans, which is the rate-controlling step of the stepwise mechanism, must show E A equal to the difference between the VGS energies of reactant and product, that we call ∆E (see Fig. 2c). The height and width of the barrier impact the absolute value of the rate, but they do not affect the Arrhenius slope in the low temperature limit of an asymmetric reaction. Indeed, a harmonic estimation of ∆E is 172 meV, which is very close to the calculated E A = 190 meV from Fig. 2a.
We further analyze how tunneling manifests itself in these reaction rates in Table I. A standard procedure to estimate the impact of tunneling is to compare k inst with the rate predicted by the Eyring TST (k TST ), since the latter neglects tunneling but includes ZPE. The tunnelling enhancement factor is κ tun = k inst /k TST [52], which is reported in Tab. I. Surprisingly, because these factors are close to 1, they would seem to indicate that tunnelling plays a minor role. To understand this observation, we computed the kinetic isotope effect (KIE), defined as k H /k D , where k D was obtained from calculations where the inner-cage hydrogen atoms were replaced by deuterium. If tunnelling would be a minor effect, the only difference in these rates should be ZPE, and since ZPE is captured by both k inst and k TST , the KIE of both should be similar. However, as shown in Tab. I, these numbers are different. TST overestimates the KIE in this particular case, for reasons outlined in the SI, Section VII. We thus conclude that κ tun ≈ 1 because of the following. On Cu(110) the ∆E (including harmonic ZPE) between reactants and products, which is a good estimate for E A in k inst as discussed above, happens to be similar to the energy difference between the ZPE corrected reactant and transition state, which defines E A for k TST (see Table II). Additionally, because close to T c both rates are comparable, if E A are similar, the prefactors also must be. This observation explains why TST fared reasonably well in these systems in the past, even without including the relevant physics of tunneling [44]. Thus, we propose that KIE inst /KIE TST can be an alternative measure of tunneling contributions to hydrogen transfer reactions.
The investigation of the dynamical surface effects on the DHT required RPI calculations where we fixed the surface atoms at the reactant position. The rates obtained as a result of this constrained optimization lack all contributions from fluctuations of the surface degreesof-freedom. We call the ratio between the rates with and without those constraints the "surface fluctuations enhancement" (SFE). The SFE for k inst on Cu(110) are reported in Table I and can adopt surprisingly large values. The results show that the dynamical surface effects act on the opposite direction of the static ones, increasing the tunneling rates up to two orders of magnitude. Interestingly, the SFE become larger at lower temperatures because the contribution of heavy atoms to tunneling increases with decreasing temperature.
We then compare the Cu(110) with the Ag(110) substrate, a surface with a weaker static interaction. On Ag(110), porphycene is less buckled upon adsorption (see Fig. 1b) and the trans conformer lies only 20 meV above cis, but T c is also 264 K. Accordingly, it was observed in experiments that the DHT rates are substantially faster on Ag(110) than on Cu(110) [4]. Unlike Cu(110), the measured rates show two distinct regimes [4], reproduced in the inset of Fig. 2b. Above ∼ 10 K there is an Arrhenius behaviour, while below ∼ 10 K the rate shows al-  Fig. 2 and text. The Tt interval is given by considering calculated or experimental references for the Ag(110) case (see SI). The surface reconstruction of Au(110) [53] was ignored for the sake of comparison. most no temperature dependence. In Fig. 2b we show the calculated k inst for the stepwise and the concerted mechanisms of the cis→cis DHT of porphycene at Ag(110) (see calculation details in the SI). We obtain an E A of 45 meV for the stepwise mechanism, which compares reasonably well with the experimental E A , which we estimate to be 12 ± 3 meV. The calculated harmonic value of ∆E is again close to E A for this reaction. At 10 K the x inst for the concerted mechanism starts at the reactant minimum, indicating that tunneling takes place from the reactant VGS. As such, the rates for the concerted mechanism (which is symmetric) do not change with temperature below this point and it becomes dominant below 8.5 K. Hence, we can explain the two regimes observed in experiment by the change in the DHT mechanism. The lack of quantitative agreement between the calculated and measured, transition temperature and rate of the concerted mechanism comes most likely from the potential energy surface used in the calculations, but some dependence of the measured rates on the STM tip cannot be fully discarded. Finally, as shown in Table I, the SFE are smaller here than they were for Cu(110), accounting for a factor 4 increase of the rate at 75 K. This is consistent with the weaker adsorbate-surface interaction strength.
Building up on these considerations, the dependence of these DHT rates can be schematically understood as shown in Fig. 2d and predicted for other metallic fcc [110] surfaces. One needs to compute the ZPE-corrected energy difference between the cis and trans tautomers ∆E, the ZPE-corrected energy barrier for the stepwise E * step and the T c of the stepwise reaction. At high temperatures the reaction behaves classically and proceeds by hopping over the lower barrier, which is normally the stepwise one, yielding a slope of ≈ E * step , labeled 1 in Fig. 2d. Considerably below T c , the low temperature limit of the stepwise tunneling reaction is achieved and a slope of ≈ ∆E should be observed (labeled 2). Finally, below a transition temperature T t , the concerted mechanism becomes dominant and the rate becomes independent on temperature (labeled 3).
Using Ag(110) as a reference and building a 1D potential model for which it is necessary to calculate the barrier for the concerted mechanism E * con (see model in the SI), T t can be estimated for other surfaces. In Table II the calculated values for different fcc [110] surfaces are reported, together with the corresponding estimation of T t . All T c values are similar and close to 300 K, showing the importance of tunneling at considerably high temperatures. While the estimated T t represent temperatures that can be achieved in different experiments especially for the stronger interacting surfaces, the resulting rates in Cu(110), Ni(110), and Pd(110) would be smaller than 10 −10 Hz, which lies far beyond the STM detection limit.
To conclude, we have shown how surface interactions can impact tunneling within a prototype molecular switch based on porphycene molecules adsorbed on metallic surfaces. This study was able to show the inadequacy of dimensionality-reduction schemes on these problems, the counter-intuitive origin of different temperature-dependencies of the rates in the deeptunneling regime, and the effects of surface interaction on the intramolecular hydrogen tunneling. Even though full-dimensional calculations are required to get quantitative results and understand the underlying processes, we could propose a simple estimator to predict the DHT temperature dependence on different metallic surfaces.
The methodology we presented can be straightforwardly applied to other molecules on surfaces where the calculation of internal hydrogen transfer rates are sought. Limitations in the RPI approximation may arise when several local minima of the adsorbate with similar energies are present [43].
The well-defined system addressed in this work allowed us to disentangle and quantify static and dynamic effects of the environment (in this case the metallic surface) on quantum hydrogen dynamics. We showed that dynamical effects of the environment can promote hydrogen tunneling. Such a quantification is normally not straightforward in condensed phase or biological systems. In this sense, this work shows how single-crystal substrates can be an ideal playground where cutting-edge theory and experiment can meet to provide a deeper understanding of quantum dynamics in fluctuating environments. These findings will help to address hydrogen dynamics in biology [54] and in functional materials [55], as well as guide the design and interpretation of future experiments.  3.501Å. See caption in Tab. S1.

II. REACTION AND BARRIER ENERGIES
In Tab. S8, the cis→ trans energy difference (∆Ẽ), the energy barrier of the stepwise mechanism (Ẽ * step ) and the energy barrier of the concerted mechanism (Ẽ * con ) for different metallic surfaces are presented. Note that, differently to the values reported in Tab

III. RING POLYMER INSTANTON CALCULATIONS
The ring polymer instanton (RPI) simulations were performed using the implementation created by us in the i-PI package [8,9]. Several enhancements in the implementation used for these calculations are available in the open-source code repository and described in Ref. [10].
The force convergence criterion for the transition state geometry and instanton optimizations was 0.005 eV/Å. The calculations need to be converged with respect to the number of discretization points of the closed closed Feynman path (CFP), also known as beads. The high-dimensional neural network reported in Ref. [11] which delivers an accurate description of the gas phase porphycene potential energy surface was used to perform such convergence tests due to its reduced computational cost. Since it presents a similar vibrational spectrum to the systems considered here, it shows a comparable convergence behaviour. In Tab. S9 and Tab. S10 the convergence of the hydrogen transfer rates for concerted and stepwise mechanisms, respectively, is presented. The reported rates for the DHT of porphycene on Cu(110) and Ag(110) in Fig. 2 in the main text were calculated with a different number of beads at each temperature, in such way to ensure that the error was always below 20% for the stepwise mechanism and within the order of magnitude for the concerted one.   As discussed in the main text, on Ag(110) at 8.5K the concerted mechanism becomes dominant over the stepwise one. The Euclidean action can be written as S = β U RP , where U RP is the potential energy of the ring polymer which represents the discretization of the closed Feynman path [12]. The exponential dependence that the rate has on S[x inst (β)] results that at low temperatures (high β values) a small error in S has a huge impact on the rate. For example, an error of only 15 meV produces an error in the rate of ≈ e 17 at 10 K, while at 100 K the error is reduced to only a factor of 5. Thus, while the results for Cu(110) are semi-quantitative, the ones corresponding to the Ag(110) surface at low temperatures should be interpreted only qualitatively.
The cis → cis reaction in the gas phase and on surfaces only admit stationary trajectories that represent second order saddle points rather than first order ones. In a similar spirit to what has been reported in Ref. [13], the steepest descent approximation along the second negative mode can be replaced by the full integration in that direction. This procedure was tested on representative 2D systems performing the required numerical integration using 200 S8 points. The results of this procedure were contrasted with taking the module of the second negative frequency and using the standard expression for the rate. While the results differ by one or two orders of magnitude, the temperature dependence is similar (same order of magnitude). Nevertheless, this effect is minor in comparison to the inaccuracy of several orders of magnitude at low temperatures described before.

VI. HEAVY ATOMS CONTRIBUTION TO THE TUNNELLING PROCESS
In Tab. S11, the tunneling path lengths, computed as the arc length along the instanton pathways, for the different atomic species at several temperatures are presented. As expected, the lower the temperature the higher is contribution from heavy atoms to the tunneling dynamics. Even though the DHT is much slower at Cu(110) than at Ag (110) where Q ‡ and Q r refer to the transition state and reactant partition functions, respectively, S9 E * to the potential energy barrier, and β = 1/k B T with k B the Boltzmann constant and T the temperature. RPI can be expressed in a similar way, see for example Eq. 1 in the main text and Eq. 8 in Ref. [15], where the terms entering the prefactor are explicitly given. In the main text the activation energy E A and prefactor refer to the effective Arrhenius slope, ∂β , and effective Arrhenius prefactor, A = k/e −βE A . In Table S12, and in order to allow a rigorous comparison between the considered rate theories, the terms "exponential factor" and "prefactor" refer, instead, to the factors that appear in Eq. 1 above and in Eq.
1 of the main text.
As discussed already in the main text, the behaviour of the DHT rate in the low temperature limit for the stepwise mechanism presents an effective Arrhenius slope with a value of approximately ∆E. In the case of the hydrogen isotopologue, ∆E happens to coincide with E * step while for the deuterium isotopologue E * step is 90 meV larger. As a consequence of these energetic relations, and a fortuitous compensation (see Tab. S12), DHT rates predicted by TST for the hydrogen isotopologue are unexpectedly similar to the ones predicted by RPI.
However, the rates for the deuterium isotopologue are largely underestimated below the cross-over temperature which results in the observed overestimation of KIE.  S12. Dimensionless exponential factor (e −S inst / and e −βẼ * step ) and prefactors (A inst and A TST ) obtained from RPI and TST simulations, respectively.Ẽ * step represents the energy barrier for the stepwise mechanism without zero-point energy corrections, because the exponential factor including these corrections is formally in A TST . Interestingly, even though both theories predict similar DHT rates in this temperature range, the corresponding factors differ significantly.

VIII. ESTIMATION OF TRANSITION TEMPERATURE FROM STEPWISE TO CONCERTED MECHANISM IN THE DEEP TUNNELING REGIME
The low temperature behaviour of the DHT rate for the stepwise mechanism can be expressed by an Arrhenius-like formula where the zero-point energy corrected reaction energy, (∆E in Tab. II.), is a good approximation of the effective activation energy, E A , and the prefactor, A step , can be approximated by the values obtained from the instanton calculations performed on the Ag(110) and Cu (110) surfaces. On the contrary, the plateau at the low temperature limit of the rate for the concerted mechanism is difficult to obtain (see Sec. V). Since this reaction is symmetric, a rough estimate of the rate can be obtained from its tunneling splitting, ∆, as A reasonable 1D approximation of the potential energy along the reactive coordinate can be written as, where V 0 is taken as the reaction energy barrier (Ẽ * con in Tab. S8), and x 0 is a measure of the barrier width that has to be determined.
In Tab. S13 and S14, the molecular distortions h bu , measured as the difference between the height of the amino and imino N atoms in the molecule, are presented for the cis tautomer  Other scaling factors based on the mass-scaled displacement between the cis adsorption conformation and a flat conformation were also tested with similar results.
Finally, the tunnelling splittings were computed using the Wentzel-Kramers-Brillouin (WKB) approximation [17], and the transition temperature T t is derived by setting k T →0 con = k T →0 step . Both values of T t , calculated from the experimentally derived parameters and from the theoretically derived parameters, are reported in Tab. II in the main text, giving the reader an assessment of the uncertainty of the estimation. The lower (higher) temperatures correspond to the calculation where the theoretical (experimental) Ag(110) rate value was