Anharmonicity in a double hydrogen transfer reaction studied in a single porphycene molecule on a Cu ( 110 ) surface

Anharmonicity plays a crucial role in hydrogen transfer reactions in hydrogen-bonding systems, which leads to a peculiar spectral line shape of the hydrogen stretching mode as well as highly complex intra/ intermolecular vibrational energy relaxation. Single-molecule study with a well-defined model is necessary to elucidate a fundamental mechanism. Recent low-temperature scanning tunnelling microscopy (STM) experiments revealed that the cis 2 cis tautomerization in a single porphycene molecule on Cu(110) at 5 K can be induced by vibrational excitation via an inelastic electron tunnelling process and the N–H(D) stretching mode couples with the tautomerization coordinate [Kumagai et al. Phys. Rev. Lett. 2013, 111, 246101]. Here we discuss a pronounced anharmonicity of the N–H stretching mode observed in the STM action spectra and the conductance spectra. Density functional theory calculations find a strong intermode coupling of the N–H stretching with an in-plane bending mode within porphycene on Cu(110).

Anharmonicity in a double hydrogen transfer reaction studied in a single porphycene molecule on a Cu(110) surface S. Liu, † a D. Baugh, † ‡ a K. Motobayashi, b X. Zhao, c S. V. Levchenko, cd S. Gawinkowski, e J. Waluk, ef L. Grill, g M. Persson h and T. Kumagai * a Anharmonicity plays a crucial role in hydrogen transfer reactions in hydrogen-bonding systems, which leads to a peculiar spectral line shape of the hydrogen stretching mode as well as highly complex intra/ intermolecular vibrational energy relaxation.Single-molecule study with a well-defined model is necessary to elucidate a fundamental mechanism.Recent low-temperature scanning tunnelling microscopy (STM) experiments revealed that the cis 2 cis tautomerization in a single porphycene molecule on Cu(110) at 5 K can be induced by vibrational excitation via an inelastic electron tunnelling process and the N-H(D) stretching mode couples with the tautomerization coordinate [Kumagai et al.Phys.Rev. Lett.2013, 111,   246101].Here we discuss a pronounced anharmonicity of the N-H stretching mode observed in the STM action spectra and the conductance spectra.Density functional theory calculations find a strong intermode coupling of the N-H stretching with an in-plane bending mode within porphycene on Cu(110).
Hydrogen (H) transfer reactions are of fundamental importance in nature. 1 The H transfer dynamics is complicated by anharmonic potential energy surfaces caused by the hydrogen bond which acts as a pathway of H atoms.The anharmonicity results in a strong intermode coupling of A-H stretching (A is an electronegative atom) with other lower frequency modes, which leads to highly complex spectral features in vibrational spectra, particularly in the A-H stretching band. 2 Although anharmonicity in H-bonding systems has been studied extensively using vibrational spectroscopy and quantum chemical calculations, 3 an accurate and quantitative description of the peculiar vibrational spectral features as well as the H-transfer dynamics remains a challenging topic.Since the anharmonic nature is hidden by inhomogeneous effects in bulk samples, it is essential to study a well-defined model and employ local spectroscopy down to the single-molecule level.Intramolecular H transfer reactions, so-called tautomerization, have served as an important model for studying a fundamental mechanism of H-transfer dynamics. 46][7][8] Recent low-temperature scanning tunnelling microscopy (STM) experiments revealed that the cis 2 cis tautomerization of a single porphycene molecule on Cu(110) can be induced by vibrational excitation via inelastic electron tunnelling and that the N-H stretching mode couples with the tautomerization coordinate. 90][11] Here we provide an in-depth analysis of single-molecule vibrational spectroscopy of STM and discuss a pronounced anharmonicity of the N-H stretching mode.

Results and discussion
Fig. 1a shows STM images before and after the cis 2 cis tautomerization of a single porphycene on Cu(110) at 5 K. 9,10 The electron-induced tautomerization of porphycene was investigated using STM action spectroscopy (STM-AS) which has emerged as a powerful method to study single-molecule reactions on surfaces. 12,13e measured STM-AS for normal (h-porphycene) and deuterated porphycene (d 2 -porphycene in which only the inner H atoms are replaced by D) 9 and spectral fitting analysis was carried out using an established method. 12,13This fitting analysis enables us to extract the vibrational properties of molecular adsorbates which are associated with the reaction.The reaction yield Y at a specific bias voltage V is defined by where R is the reaction rate, I t the tunnelling current, and e the elementary charge.At constant V, R follows a power law dependence, 14 where N is the total reaction order, i.e., the number of tunnelling electrons required to trigger the reaction.Theoretically the reaction rate is given by the vibrational generation rate (G iet ) via an inelastic electron tunnelling process. 15G iet can be expressed as a convolution of the vibrational density of states (DOS, r ph (o)) and the spectral generation rate function (G in ): G iet ðVÞ ¼ Ð 1 0 dor ph ðoÞG in ðoÞ, as derived in ref. 15.Y(V) can be expressed using double integration of the vibrational DOS as where K is the prefactor (rate constant) that is determined via the elementary process behind the reaction and n is the reaction order of vibrational excitation.If several vibrational modes are involved in STM-AS, the total reaction yield Y tot (V) becomes the summation of Y(V) from each active vibrational mode i with energy hO i , 12,13,16 Y tot ðVÞ ¼ When the reaction is induced by simultaneous excitation of several different vibrational modes, the total reaction order (N) involves contribution from each excitation.8][19][20][21] Here we use a Gaussian function r ph ðoÞ ¼ ! because it allows taking into account broadening effects including all possible contributions, i.e., finite temperatures, vibrational relaxation, instrumental resolution, 12 and anharmonic effects due to hydrogen bonding. 16The full width of half maximum (g i ) of the Gaussian vibrational DOS is given by g i ¼ 2 ffiffiffiffiffiffiffiffiffi 2ln2 p s i .We first discuss the STM-AS for the cis 2 cis tautomerization of d 2 -porphycene (Fig. 1b).The lateral tip position during the measurement is indicated in Fig. 1a, which corresponds to the side of pyrrole rings with the amine N atoms.Fig. 1c shows the current dependence of the tautomerization rate, where N changes from B2 to B1 around 260 mV.Thus, tautomerization occurs via a two-electron process and a one-electron process at the low and high voltage regimes, respectively.At higher currents, the tip-molecule distance becomes smaller and the STM tip may affect the tautomerization through the interaction with the molecule, which eventually leads to modification of the potential energy surface. 11However, the modification is expected to be rather small at the tip-molecule distances in this experiment (the tip is substantially far from the equilibrium distance of the tip-molecule potential, where the significant deformation of the potential occurs 11 ).Hence, the influence on the reaction order should be negligible.Since the tautomerization rate is very small at low bias voltages, we used a relatively high current to measure the yields (see the lower panel of Fig. 1b) in order to observe a sufficient number of tautomerization events in a realistic time-scale of the experiment.However, according to eqn ( 1) and ( 2), Y(V) depends on the current in a multi-electron process, in contrast to a one-electron process.Therefore, the measurement current affects the Y(V) curve at low voltages.The small markers in the upper panel of Fig. 1b show the measured (raw) yields (Y raw ) obtained at the used tunnelling current (I raw ) indicated in the lower panel.In order to examine the influence of the current to the Y(V) curve, we test the following normalization.The open markers in Fig. 1b are the yields (Y norm ) normalized to a specific current (I norm : 0.1, 1, and 10 nA), which are given by raw Y raw where N was obtained from a fitting analysis of the current dependence measurement (cf.Fig. 2b).These values for I norm are in a typical range of tunnelling currents in STM.As mentioned above, we do not take into account any influence of the tip-molecule interaction that may modify the potential energy surface of the tautomerization. 11Note that this normalization was not carried out in our previous report. 9As discussed later (cf.Fig. 2c), N may be varied depending on the measured current range when one-and two-electron processes are simultaneously involved.The consistency of the normalized yield was confirmed at several bias voltages where the current dependence of the rate is available (Fig. 1c).It is clear that Y(V) at different currents exhibits a vertical shift in the two-electron regime (Fig. 1b), which leads to a critical change in K i in eqn (5)  but the other vibrational factors, namely the vibrational energy and width, are not significantly influenced in the fitting analysis.For d 2 -porphycene, Y(V) exhibits a threshold at B160 mV and a second steep increase at B260 mV.A similar result was obtained at both bias polarities.The solid line in Fig. 1b is the best-fitting result to eqn (5) with the parameters listed in Table 1.The obtained vibrational energy for the second increase of Y(V) is in excellent agreement with the calculated symmetric (anti-symmetric) N-D stretching mode, n s(as) (N-D), of 281.9 (280.9)meV.However, the vibrational DOS around the threshold voltage show a large broadening factor which ranges typically from 5 to 15 meV. 13This broad vibrational DOS may be explained by the contribution from several different vibrational excitations.Additionally, the reaction is second order (i.e., twoelectron process) at this voltage range, thus multiple vibrational excitations are also involved and the reaction should comprise complex elementary processes. 18Moreover, a bunch of vibrational modes exist for porphycene below 200 meV, 9 hampering a clear assignment of the involved modes.A similar threshold voltage (B150 mV) is also observed for the tautomerization of porphycene on Cu(111). 22,23A recent theoretical study by Novko et al. proposed that the excitation of skeletal modes plays a crucial role in the reaction. 24It should also be noted that the reaction (minimum energy) path for tautomerization on Cu(110) and Cu(111) does not initially involve N-H stretching but rather a small lateral translation of porphycene along the surface. 10,11,25Therefore, the reaction coordinate should also involve low frequency modes such as the hindered translation and rotation of the molecule.
The tautomerization is observed to be enhanced by the n(N-D) excitation, which suggests that this reaction does not follow the reaction path since this vibrational energy does not appear in low frequency modes along the initial reaction path.7][28] K i in eqn (5), in principle, contains important physical parameters such as the vibrational dumping rate of the excited vibrational mode through relaxation to the reaction coordinate and/or other low frequency modes (including surface phonons) and electron-hole pair excitation. 29,30However, we were unable to extract such valuable information due to the lack of the detailed knowledge of the complex elementary process of tautomerization of porphycene on the surface.
Y(V) of h-porphycene shows a conspicuous difference from d 2 -porphycene (Fig. 2a).Although the threshold at B150 meV is similar to that of d 2 -porphycene, the second increase of Y(V), which should correspond to the n(N-H) excitation, is not as obvious as n(N-D) excitation of d 2 -porphycene.Instead, Y(V) exhibits a very moderate increase ranging from 200 to 380 mV.A similar moderate increase of Y(V) was also observed for the H-bonded O-H stretching mode in a water-hydroxyl complex on Cu(110), which is interpreted as vibrational broadening caused by an anharmonic potential. 16As shown in Fig. 1b, Y(V) of d 2 -porphycene can be reproduced with eqn (5) using two vibrational DOS, that is, excitation of skeletal modes at 181.3 meV and n(N-D) at 274.7 meV.Using a simple harmonic approximation, the n(N-H) energy can be estimated to be 375 meV from the n(N-D) energy, 31 which is in excellent agreement with the simulated n(N-H) at 376 mV by DFT. 9 However, if the STM-AS of h-porphycene is simulated with eqn (5) using two vibrational DOS, namely skeletal mode around the threshold voltage and n(N-H) at 375 meV with a comparable broadening factor obtained for d 2 -porphycene, the resulting curve largely deviates from the experimental result (see the red dashed line in Fig. 2a).
A remarkable difference between h-and d 2 -porphycene is also observed in the voltage dependence of N (Fig. 2b).The transition from a two-to a one-electron process is relatively sharp for d 2 -porphycene, whereas it is much more gradual for h-porphycene 32 and the transition spreads from 220 to 270 mV.This gradual transition suggests that the vibrational DOS of the skeletal mode and the N-H stretching are simultaneously excited in this voltage range, resulting in the non-integer N. In this situation, the total reaction rate is given by 19 where k 1 [s À1 A À1 ] and k 2 [s À1 A À2 ] are the rate constants for the one-and two-electron processes, respectively.In Fig. 2c, the current dependence of the tautomerization rate obtained at V sample = À250 mV is fitted by eqn ( 6), giving k 1 = 4.99(AE0.73)Â 10 6 s À1 A À1 and k 2 = 3.02(AE0.30)Â 10 17 s À1 A À2 .At low currents the one-electron process dominates the total rate, whereas the two-electron process prevails at high currents.Therefore, N in eqn ( 2) is affected by the measured current range when different processes compete.We also note that the second term in eqn (6) should be more complex because the two-electron process may involve various elementary processes. 19,33As discussed above, the vibrational DOS around the threshold voltage involve the contribution from multiple vibrational excitations.Thus, an in-depth analysis of k 2 requires precise information about contributions from the individual modes, which is not available in the present case.A detailed theoretical analysis for a complex process including one-and two-electron processes has been applied only to a simple reaction such as rotation of an acetylene molecule on Cu(100). 34 remarkable difference between h-and d 2 -porphycene was also observed in the dI/dV spectrum, 9 where vibrational signals appear as a characteristic peak or dip (Fig. 3a).The peak and dip at 177 and 284 mV observed for d 2 -porphycene nicely match the skeletal modes and the n(N-D) excitation, respectively.These values are also consistent with the fitting parameters of Y(V) (Table 1).However, the dI/dV spectrum of h-porphycene exhibits rather complex features between 300 and 400 mV consisting of a broad peak and a dip.Furthermore, it is found that the dip includes two components in the high-resolution dI/dV spectrum (see the inset of Fig. 3a).These spectral features suggest the existence of several vibrational excitations around n(N-H).We find that Y(V) of h-porphycene can be reproduced with eqn (5) using at least two vibrational DOS in the n(N-H) a These values are restricted around the peak and dip position observed in the dI/dV spectrum of h-porphycene in Fig. 3a.
regime where the peak and dip are observed in the dI/dV spectrum.The fitting result is indicated by the solid line in Fig. 2a and the fitting parameters are listed in Table 1.
The n(N-H) energy found in the Y(V) and dI/dV spectra appears to be considerably red-shifted as compared to the N-H stretching mode, for example, in NH 3 (417 meV 35 ), freebase porphyrin in rare-gas matrices (412 meV 36 ), and phthalocyanine on Ag(111) (406 meV 37 ), indicating the presence of relatively strong H bonds in porphycene on Cu(110).][44] Since no fundamental mode exists between 200 and 376 meV for h-porphycene on Cu(110), 9 the peak at 330 mV in the dI/dV spectrum observed for h-porphycene should be assigned to overtone/combination excitation of lower vibrational modes.The absence of a peak (observed for h-porphycene) around the n(N-D) energy region of d 2 -porphycene implies that the peak at 330 mV for h-porphycene appears in a similar mechanism to other H-bonding systems (as discussed above) and the coupling may be enhanced through a Fermi resonance between n(N-H) and the overtone/combination mode(s), most likely N-H bending.Additionally, it is found that the second dip at 385 meV observed for h-porphycene is absent in d 12 -porphycene in which all the peripheral H atoms are replaced by D atoms but without any substitution of the inner H atoms (see the chemical structure in Fig. 3b).This indicates that the first and second dips cannot be assigned to symmetric and anti-symmetric combinations of the fundamental n(N-H).The possibility of C-H stretching modes is also excluded because the dip is also absent in d 2 -porphycene.Therefore, it should also be assigned to the different overtone/ combination modes from the one at 330 mV.
The spectral anomalies of n(N-H) of porphycene on Cu(110) indicate a pronounced effect of the anharmonicity causing a strong coupling with a lower vibrational mode(s) that can modulate the H-bonding geometry within the molecular cavity.Such a coupling mechanism and vibrational line broadening of n(N-H) were examined with quantum mechanical calculations for porphycene in the gas phase. 5Here we have investigated one anharmonic coupling term on the surface by calculating the variation of the n(N-H) and n(N-D) energies as a function of the displacement along a lower frequency mode (Fig. 4c).This mode causes the largest displacement of the N-HÁ Á ÁN geometry and is essentially an N-H bending mode (Fig. 4b).A considerable variation occurs in the n(N-H) energy (Fig. 4c), indicating a strong anharmonic coupling of n(N-H) with the bending mode.However, this variation is almost negligible for n(N-D) because the anharmonicity of the potential has lesser impact on n(N-D) than on n(N-H) because of the smaller zero-point motion of D than H. Finally, note that our study only indicates that there is a strong anharmonic coupling term for h-porphycene compared to d 2 -porphycene but is not intended to describe the anharmonic vibrational line shape and the Fermi resonance.A more complete description would require the calculation including more anharmonic terms in the coupling between the stretching mode and the bending mode as well as nuclear quantum effects.

Conclusions
We presented a detailed analysis of STM-AS for the cis 2 cis tautomerization of porphycene isotopologues on a Cu(110) surface at 5 K and discussed a pronounced anharmonicity of the N-H stretching mode.The spectral fitting analysis of STM-AS in combination with the reaction order measurements revealed the vibrational energies and broadening factors of the involved modes.The tautomerization is induced through multiple excitations of skeletal vibrational modes of porphycene and the N-H(D) stretching largely enhances the reaction.These vibrational excitations were also identified in the conductance spectroscopy, manifested as a peak or dip around the vibrational energy.We found that the N-H stretching mode shows not only a significant broadening, but also an emergence of additional vibrational features which indicate a pronounced anharmonicity and strong coupling with lower frequency modes.In contrast, these spectral anomalies were absent in the N-D stretching mode.Theoretical simulations demonstrated a strong intermode coupling of the N-H stretching with the bending mode but not for the N-D stretching mode.These results indicate a significantly different contribution of the anharmonicity between the N-H and N-D stretching modes.Our approach paves the way for studying anharmonicity in H-transfer reactions on solid surfaces at the single-molecule level.

Experiments
All experiments were performed in an ultra-high vacuum chamber (base pressure of 10 À10 mbar) equipped with a lowtemperature STM (modified Omicron instrument with Nanonis Electronics).STM measurements were conducted at 5 K and the bias voltage (V bias ) was applied to the tip (V tip ) or sample (V sample ).The Cu(110) surface was cleaned by repeated cycles of argon ion sputtering followed by annealing to 700 K. STM tips were made from a tungsten or gold wire and then optimized in situ by applying a voltage pulse and controlled indentation of the tip into the surface.Porphycene molecules were deposited onto the surface at room temperature from a Knudsen cell (at a temperature between 450 and 500 K) and the sample was then transferred to the STM at 5 K.

Computational methods
DFT calculations were performed using the Fritz-Haber Institute ab initio molecular simulations package (FHI-aims) 45 with the Perdew-Burke-Ernzerhof (PBE) functional 46 and a Monkhorst-Pack 47 k-grid of 3 Â 3 Â 1 k-points in the supercell.Long-range van der Waals forces were accounted for by the Tkatchenko-Scheffler scheme. 48An all-electron and default tight basis set was used.The Cu(110) surface was represented by a four-layer slab with a 4 Â 6 surface unit cell and a 420 Å vacuum region.The two topmost layers were relaxed and the residual Cu atoms were fixed at the bulk geometry (lattice constant of 2.568 Å).The vibrations of the adsorbed molecule were calculated on a rigid substrate lattice (the Cu atoms were fixed) by diagonalising the dynamical matrix which was obtained via finite differences of the calculated forces at symmetric ionic displacements of 0.02 Å.

Fig. 1
Fig. 1 (a) STM images of a single porphycene molecule on Cu(110) at 5 K (V tip = À100 mV, I t = 10 nA, size: 1.49 Â 1.42 nm 2 ). 10 The white star indicates the tip position during the STM-AS and conductance measurements.The white grid lines represent the surface lattice of Cu(110).(b) STM-AS for d 2 -porphycene on Cu(110) measured at 5 K.The small closed circles and squares in the upper panel indicate the raw data obtained at the tunnelling current shown in the lower panel.The open markers in the upper panel are the calculated Y(V) normalized to a constant current indicated in the figure.The solid curve represents the best-fitting result of Y(V) for I t = 1 nA to eqn (5) with the parameters listed in Table 1.(c) Current dependence of the tautomerization rates (R H-L ) measured at various bias voltages.The slopes (N) are determined to be 1.83(AE0.11),2.12(AE0.13),1.95(AE0.08),2.12(AE0.13),1.03(AE0.07),0.88(AE0.06)and 0.71(AE0.05)at 200, 240, 250, 265, 280, and 350 mV, respectively.

Fig. 2
Fig. 2 (a) Y(V) for h-and d 2 -porphycene on Cu(110) measured at 5 K.The yields in the N 4 1 regime are normalized to a constant current (1 nA).The dashed red line represents the curve obtained from eqn (5) using the parameters: (h O, s ph ) = (181 meV, 23.2 meV) and (375 meV, 6.8 meV).The solid red line represents the best-fitting result with the parameters listed in Table 1.The multiple data points at several bias voltages represent the yields obtained either at a different current or under different tip conditions.(b) Voltage dependence of the total reaction order (N) for h-and d 2 -porphycene.The experimental data (makers) are fitted to a step function convoluted with a Gaussian profile.The broadening factors (FWHM) of the Gaussian function are B42 and B12 meV for h-and d 2 -porphycene, respectively.(c) Current dependence of the tautomerization rate measured for h-porphycene at V sample = À250 mV.The black solid line is the best fit using eqn (6), and the dashed lines indicate the contribution from oneand two-electron processes.The inset graph plots the logarithm of the rate and current.

Fig. 3
Fig. 3 (a) dI/dV spectra for h-, d 2 -and d 12 -porphycene.The spectra over the molecule are subtracted with a background spectrum measured over the Cu(110) surface.The tip was fixed at the position indicated in Fig. 1a with the gap conditions of V sample = 100 mV and I t = 20 nA.The spectra were recorded using a lock-in amplifier with a modulation voltage of 12 mV at 710 Hz frequency.The upper inset shows the high-resolution dI/dV spectrum for h-porphycene recorded with a modulation voltage of 3 mV.(b) Chemical structures of d 12 -porphycene.

Fig. 4
Fig. 4 (a) Schematic of anti-symmetric and symmetric n(N-H).(b) Calculated bending mode at 193.5 (186.6)meV for h-(d 2 -)porphycene on Cu(110).(c) Variation of the n(N-H) and n(N-D) energies as a function of the displacement along the bending mode.The horizontal axis represents the normal coordinate of the bending mode, which is normalized by its root-mean-square vibrational amplitude of the ground state.

Table 1
Fitting parameters of Y(V) for h-and d 2 -porphycene