Yi-Han Cheng(程奕涵), Yu-Cheng Zhu(朱禹丞), Xin-Zheng Li(李新征),2,3,4, and Wei Fang(方為)

1State Key Laboratory for Artificial Microstructure and Mesoscopic Physics,Frontier Science Center for Nano-optoelectronics and School ofPhysics,Peking University,Beijing 100871,China

2Interdisciplinary Institute of Light-Element Quantum Materials and Research Center for Light-Element Advanced Materials,Peking University,Beijing 100871,China

3Peking University Yangtze Delta Institute of Optoelectronics,Nantong 226010,China

4Collaborative Innovation Center of Quantum Matter,Peking University,Beijing 100871,China

5Department of Chemistry,Fudan University,Shanghai 200438,China

6State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical Computational Chemistry,Dalian Institute of Chemical Physics,Chinese Academy of Sciences,Dalian 116023,China

Keywords: quantum tunneling,proton transfer,hydrogen bonding

1. Introduction

The abundance of hydrogen(H)makes H-bonding ubiquitously present in a wide range of systems, such as in liquid water, ice, proteins, just to name a few.[1–10]H-bonding networks are among the most important and interesting structures in nature, for example, they can give rise to a variety of condensed water phases; they hold proteins in their functioning structure, as well as creating proton conduction chains. The small mass of H means that many properties of the H-bonded systems are quantum mechanical in nature.[11–14]The associated nuclear quantum effects(NQEs),which refer to the differences between properties of a realistic system when the nuclei are described as quantum or classical particles,therefore,lie at the heart of accurate simulations of the physical and chemical properties of such systems. These NQEs can be reflected by many phenomena,e.g.,quantum delocalization of H,quantum tunneling,and zero-point energy(ZPE)effects.[15–19]

Proton transfer along H-bonds are a class of processes that possess rich quantum behavior. While it is fair to say that in processes involving a single PT, the microscopic details have been relatively well-understood,[20–26]when multiple protons are involved,however,the underlying mechanisms become more complicated and atomic-scale understanding is still lacking. One prototype example for the existence of multiple PT processes resides in systems where multiple H-bonds are arranged in a cyclic manner.In such systems,it can happen that each proton moves separately meaning that the PT process occurs in a“stepwise”manner,where multiple metastable intermediate states (ISs) typically exist, resulting in multiple transition states (TSs). It can also happen that the protons transfer collectively,commonly referred to as the“concerted”mechanism,and the process has only a single TS.For systems with multiple meta-stable ISs,that the potential energy barrier for stepwise transfer is usually smaller, so that some people even believed that the stepwise PT pathway would always be dominate.[27]Whether the concerted transfer mechanism can be important has been under debate.[28]However, more recently, concerted transfer has been proposed to exist as well in several systems,such as systems with short H-bonds,which may result in low reaction barriers for concerted PT.For systems without meta-stable ISs,typically,stepwise transfer creates charge defects,leading to an unfavorable potential energy barrier,which hinders this mechanism,while concerted transfer circumvents the associated energy penalty.[19,29]Therefore,from comparing the barrier heights of the two pathways on the Born–Oppenheimer potential energy surface (PES), deciphering which mechanism is favored in PT processes should be straightforward. However,in-depth analysis of the experimental findings on such systems has presented a task which is much more challenging than this expectation.

A very important and practical measure to partially quantify and to probe NQEs on PT processes is isotope substitution,e.g.,replacing H by deuterium(D)or tritium(T),or replacing16O with18O,etc. Taking theoretical/experimental descriptions of the chemical reaction rates as an example,the kinetic isotope effect (KIE) is probably one of the most often used quantities when one wants to decipher the importance of NQEs in an experimental observation. These KIEs originate from both quantum tunneling and ZPE.[30,31]They are the ratio of the reaction rate of two isotopologues or isotopomers of a reaction (in most cases H in the molecules replace by D), with the rate of the lighter isotopologue divided by the rate of the heavier one. Large KIE is the main experimental and theoretical indication that atom tunneling is happening. It increases with decreasing temperature(T).[32–34]This is especially true in the deep tunneling regime,and a prominent example of this resides on H/D diffusion on metal surfaces, where quantum tunneling of H/D plays a central role.[20–26]When multiple protons are involved,the experimental KIEs are intriguing and difficult to decipher,and an atomic-scale understanding of the KIEs is still lacking.[19,35–39]

Apart from the chemical reaction rates and the KIEs,tunneling splitting may also serve as an alternative experimental observable to characterize the PT processes. Unlike rate,tunneling splitting is aT-independent property. In a symmetric double-well system, tunneling from one well to another means that the nuclear wavefunction distributes in both wells.This will mix the original degenerate states and result in a splitting of them. The mixing between the vibrational states of the degenerate potential wells leads to a splitting pattern which can typically be observed by high-resolution spectroscopic methods.[40]For systems in which the two degenerate wells are connected by a multiple PT process,it is commonly believed that ground state tunneling splittings, if measurable,serve as strong evidence for the concerted PT mechanism.[41]

Experimentally, development of state-of-the-art techniques means that one can perform isotope substitution to measure the KIEs using nuclear magnetic resonance (NMR)in molecular systems,[42–45]high-resolution scanning tunneling microscope (STM) on surfaces,[46–52]helium scattering on surfaces,[26,53]or sum-frequency generation (SFG) spectroscopy for interfaces.[54–56]Alternatively, the influence of isotope substitution on other physical properties such as the structural ones can also be measured by neutron scattering experiments[7,57]or x-ray diffraction in solids,[58,59]providing useful supplementary information to the experimental data on KIEs. Besides these, other experimental techniques such as deep inelastic neutron scattering(DINS)are also developed,which is capable of probing the proton and oxygen momentum distributions and quantum kinetic energies,intriguing new observations of the NQEs.[60,61]These experiment data provide valuable benchmark for the calibration of the theoretical methods and models.[62–70]

From the theoretical perspective, a variety of schemes can be employed to investigate such NQEs of H-bonds or KIEs of the PT processes. In the 1980s, pioneering pathintegral simulations of materials were performed with pathintegral molecular dynamics (PIMD),[71–73]where the interatomic interactions were described with empirical potentials.Later,with applications ofab initioelectronic structure calculations, especially density-functional theory (DFT), in chemical reactions,ab initio-based path-integral approaches became possible.[74–76]In recent years, advances in the pathintegral methodology have made it more feasible to study NQEs in more complicated systems[14,77–83]Concerning the mechanism of multiple PT in cyclic H-bonded systems,early studies show that in most cases, reactions involving multiple PT are not synchronous,[27,62,84]stepwise mechanism is dominant. This is especially true at highTs, where multiple PT in H-bonded networks occurs once a time via thermal activation.As a result,a charge defect moves in the H-bond network.While at lowTs hopping over the energy barriers is avoided,quantum-mechanical tunneling plays an important role. This can result in a fact that the mechanisms of multiple PT change at differentTs, and a dominance of concerted PT at lowTs is possible in some systems.[85]So far, it is fair to say that the intricate interplay between the quantum tunneling nuclei and the high-dimensional nature of these processes makes pure analysis of the Born–Oppenheimer PES inadequate to unravel the microscopic details of the PT processes in such cyclic Hbonded systems. A complete review on the experimental and theoretical studies of this topic,with feasible insights provided by new simulations,is highly desired.

In this manuscript, we present such a review, supplemented with new theoretical results. This review is organized as the following: First, we introduce theoretical methods to distinguish concerted and stepwise mechanisms. Then, we review previous experiments and simulations of multiple PT.Two kinds of systems are reviewed in this article: (i)systems with meta-stable ISs, which feature both concerted and stepwise mechanisms;(ii)systems without meta-stable ISs,which feature a concerted tunneling mechanism. According to the literature, it is found that for some systems, the discrepancy between experiment and simulation could not be neglected.Finally,we discuss some possible factors causing the discrepancy and hope new experiments/theory could settle this problem.

2. Theoretical tools for understanding PT

In this section,we introduce some of the theoretical methods for estimating tunneling rates, KIEs, and tunneling splittings. We focus on path-integral based methods, while thorough reviews for other methods can be found elsewhere.[86]

2.1. Rate theories

Classical transition-state theory (TST)[87]has been tremendously successful in predicting the rate of chemical reactions, as well as for providing theoretical insights into the underlying mechanism. In classical (harmonic) TST, one searches for a minimum energy transition pathway from the reactant to the product,and the highest energy configuration on the optimal path is called a transition state(TS).The classical TST rate at a given temperature is predominately determined by the energy(with respect to the reactant)of the TS,and the temperature dependence of the classical rate for a given reaction follows the Arrhenius behavior.Classical TST,along with methods for obtaining TSs,are among the most widely applied theoretical tools in computational chemistry.

Since classical TST does not take NQEs into account, it fails to describe systems where NQEs are important. Many quantum corrections have been developed attempting to“dress up” classical TST as a quantum rate theory.[88]For example,(harmonic)ZPE corrections,parabolic-top tunneling corrections, one-dimensional tunneling corrections based on the minimum energy pathway(MEP),etc. These corrections have also been widely applied and are successful at improving the predictive power of classical TST. However, no matter how“quantum” they seem, they still rely on the core of classical TST,which is the MEP.Tunneling pathway can deviate from the classical MEP, an effect know as corner-cutting or curvature effect. These quantum corrections fail to model complicated tunneling scenarios where corner-cutting effects are drastic, which can be problematic for studying multiple PT problems.

Feynman’s path-integral framework provides an alternative form of quantum mechanics.[89]Instead of using wavefunctions, it utilizes the classical concept of paths, making it more applicable to molecular simulations. A variety of path-integral based methods have been developed for simulation of NQEs, for example, path-integral molecular dynamics samples the quantum partition function with Newtonlike dynamics, utilizing the isomorphism between the quantum partition function and the classical partition function of a ring polymer. Comprehensive reviews of path-integral theory,its implementation into computer codes,and applications in chemistry, physics, and materials science can be found in the literature.[13,72,74,90]Considerable progress has also been made with developing quantum or semiclassical rate theories based on the path-integral framework, providing alternative routes to incorporating tunneling,other than as corrections to the classical TST.For example,the well-known ring-polymer molecular dynamics (RPMD) rate theory[91–93]has shown great performance on gas-phase reactions of small molecules compared to exact quantum mechanical benchmarks,[94,95]while being applicable to condensed phase systems.[96]

Thanks to Feynman’s path-integral, a semiclassical rate theory that formally resembles the classical (harmonic) TST has been developed, known as the ring-polymer instanton theory.[97–101]It’s an alternative formulation of the instanton theory proposed by Miller,[102]made practical for molecular simulations, and is closely related to RPMD rate theory.[98]The instanton rate is computed based on finding an optimal tunneling pathway(instanton),and is formally given as a prefactor(Ainst)times an exponential term,

Comparing with the classical TST rate, the exponential part in the instanton rate is determined by the Euclidean action of the instanton (Sinst) instead of the potential energy barrier.Thus, non-Arrhenius behavior manifests as the instanton action changes with temperature. Ring-polymer instanton theory offers a good balance between accuracy and computational cost,making it a well-rounded tool for studying tunneling processes, especially when combined withab initioelectronic structure.

2.2. KIEs and the Bell–Limbach model

KIEs can be directly simulated via direct rate calculations with isotope substitution. Alternatively, theories aimed at qualitatively predicting KIEs have also been proposed.The Bell–Limbach model is a relatively simple and easy-toapply model that links qualitative behaviors of KIEs with whether the underlying multiple PT mechanism is stepwise or concerted.[103,104]Here we summarize the idea and the conclusions of the model in simple terms(rather than rigorously)on a double proton transfer(DPT)process.

For a stepwise DPT process with a metastable IS, the Bell–Limbach model assumes that thermal equilibrium is reached at the metastable IS well. If one proton is replaced by a deuteron,the step involving the transfer of the deuteron becomes the rate limiting step in the total DPT process,therefore the total rate drops significantly. Next,if the second proton is also replaced by a deuteron,the rate is roughly halved,as the process takes to identical rate limiting steps to complete.Thus,for stepwise mechanism with IS,the KIE trend is

For concerted DPT,in the semiclassical picture, the tunneling rate is predominately determined by the abbreviated action of the tunneling path,

whereμis the effective mass of the system,Vis the potential energy along the path andE0is the reactant energy. Upon isotope substitution,the effective mass of the system changes and therefore the abbreviated action satisfies the following relation,WHH+WDD≈2WHD. Hence the KIEs qualitatively satisfy,

The above is a simple hand-waving introduction of the Bell–Limbach model, a more detailed and more rigorous analysis can be found in Ref.[103].

The Bell–Limbach model have also been extended to describing KIEs in tunneling processes involving more than two protons.[105]Here we skip the derivation and just present some of the conclusions relevant to the systems we review. For the quadruple proton transfer system, the relation between transfer rates with four transferring H atoms(k4H),three H and one D(k3H1D),and four D atoms(k4D)qualitatively satisfy

if the underlying mechanism is concerted.

2.3. Tunneling splittings

Similar to rates and KIEs, tunneling splitting levels can also be computed with methods based on the path-integral framework. In a symmetric double-well system, overlap of the wavefunction of the two wells by tunneling will result in splittings of levels in energy representation. The splitting size of ground-state is determined by the tunneling probability and reveals information on the tunneling processes that occurs in the system. Tunneling splitting can be computed with either PIMD based methods[106]or instanton based methods,[107]and we briefly review the latter in this section.

Tunneling splitting (Δ) can be expressed using quantum partition functions in theT →0 limit,[107]

whereQ(β) and 2Q0(β) are the partition functions of the double-well system and two unconnected single wells respectively.E0±Δ/2 are the energy levels split by tunneling.Qcan be expressed using kink instantons,i.e.,linear polymer instead of ring polymer instantons. The final expression for the ground-state tunneling splitting of a double-well system,

Skinkis the Euclidean action of the kink described by a linear polymer, andΦis a term related to the fluctuations around the instanton path. Both quantities areT-independent. For detailed derivation of this method,please refer to Ref.[107].

3. DPT in systems with meta-stable IS

3.1. Porphyrin and porphycene

Porphyrin is an important molecule with a wide range of applications,e.g., it has been used in designing advanced materials such as organic metals,molecular wires,and photosensitizers.[108–111]It has also been used in medicine for the treatment of cancer and dermatological diseases.[112]Free base porphyrin has two hydrogen atoms in the inner part of the skeleton. The inner-cage hydrogen atoms can adopt different configurations,giving rise to two stable tautomers,namelytransandcis(as depicted in Fig. 1), with thetranstautomer being the most stable. The two inner hydrogen atoms can hop between the different N sites, featuring a signature intramolecular DPT that can result in the switching between the tautomers. The DPT reaction connects the two degeneratetranstautomeric states can happen through a stepwise or concerted mechanism(Fig.2). This process is of vital importance in photosynthesis and metal coordination chemistry.[113]It has also made prophyrin-like molecules ideal platforms for probing multiple PT processes both experimentally and theoretically. Illustrations of the two competing DPT mechanisms(stepwise and concerted PT) in prophyrin-like molecules are shown in Fig.2. In the former case,the protons are transferred sequentially through an metastable IS (cistautomer) and the reaction pathway involves a first-order saddle point. While in the latter case,two protons are transferred in a concerted fashion, passing through a second-order saddle point. The concerted pathway in prophyrin-like molecules is prohibited classically,as the second-order saddle point is not a TS.

For porphyrin,tautomerization rates and KIEs have been measured by NMR experiments in the solvated state and in the solid state[67,115]over a range of temperatures(from 100 K to 368 K).Rate constants of the tautomerization process could be measured on the second time-scale by magnetization transfer techniques and on the millisecond time-scale using lineshape analysis.[116]Spontaneous tautomerization seizes to be measurable at low temperature. The experimental rates also show obvious non-Arrhenius behavior,[115]indicating that quantum tunneling plays an important role. The measured HH/HD/DD KIEs are consistent with Bell–Limbach model for the stepwise DPT mechanism.[115]Therefore, based on the above observations, previous studies have concluded that the DPT mechanism in porphyrin is stepwise. The metastablecistautomer, however, was not directly observed in experiments.It has been argued that the life-time of thecistautomer is too short to be measured by NMR. Theoretical results support the experimental finding that stepwise DPT in porphyrin is favorable.[117,118]DFT calculations(at B3LYP/TZ2P level)show that for porphyrin,the barriers for DPT is 16.2 kcal/mol for stepwise DPT, while the barrier for concerted PT is even higher (24.4 kcal/mol).[114]They reasoned that the high barriers are the result of the relatively weak intramolecular hydrogen bonds (as indicated by the relatively long distance of～2.9 ?A between adjacent nitrogen atoms). Therefore, based on the previous experimental and theoretical studies,the DPT mechanism in porphyrin is quite well understood, which occurs via the stepwise mechanism.

Fig.1. Structures of porphyrin. Left: trans porphyrin. Right: cis porphyrin.Reprinted with permission from Ref.[114].Copyright 1997 Springer Nature.

In 1986, porphycene has been synthesized, it is the first synthesized structural isomer of free base porphyrin.[3,119]The main difference from porphyrin is that porphycene is “rectangular” instead of “squaric”, with an NN distance of 2.6 ?A along the short side, shorter than the NN distance of 2.9 ?A in porphyrin. In drastic contrast to porphyrin, DPT phenomena in porphycene is much more complicated,despite that the two molecules are similar. Experimental characterization of the DPT in prophycene has been a long standing challenge,and in the past decades, a variety of experiments have been extensively conducted to investigate this system.[120–126]Yet,the previous experimental studies have not drawn a firm conclusion on whether the underlying mechanism is concerted or stepwise. Experiments at cryogenic temperatures on gas phase prophycene measured a tunneling splitting of 4.4 cm-1in the ground vibrational state,[127–129]indicating thattrans–transtautomerization can occur via concerted tunneling,as it is generally believed that stepwise tunneling does not result in tunneling splitting. Tautomerization rates for prophycene were measured in NMR experiments in the solid phase,which are significantly higher than that in porphyrin.[120,121]These early NMR studies could not determine whether the thermal DPT mechanism is stepwise or concerted,and there seems to be a discrepancy[130]between the observed tunneling splitting and the observed rate constants. Further experimental investigations revealed that the environment, especially the solid environment,have an impact on the PT process as well as the proton ordering in prophycene,[123]which was not the case for porphyrin. As the temperature decreases, a phase transition from a proton disordered phase I to an ordered phase II occurs at 225 K, and the DPT mechanisms as well as rates differ for the two phases. Based on the KIEs, the authors believed that the DPT in phase II occurs via the stepwise mechanism.[123]However,experimental studies with a different technique found that the DPT in porphycene is not sensitive to the environment, and that the concerted mechanism is favored due to tunneling.[126,131]

Fig.2. (a) Two-dimensional illustration of the stepwise tautomerization pathway in porphyrin. (b) Two-dimensional illustration of the concerted tautomerization pathway in porphycene. Adapted with permission from Ref.[131]. Copyright 2007 American Chemical Society.

From the theoretical perspective, various simulations have also been performed, aiming to understand DPT in prophycene and to explain the different experimental observations.[130,132–136]The presence of strong intramolecular H-bonds in prophycene (along the short side of the rectangle formed by the four N sites)results in relatively low reaction barriers for both the concerted and stepwise PT mechanism, which are 7.6 kcal/mol and 4.9 kcal/mol respectively(computed at the B3LYP/TZ2P level).[132]Since the stepwise pathway has the lower barrier,this result(from an early work)suggests that the stepwise pathway is favored. In a later study,Smedarchinaet al.[137]constructed a multidimensional Hamiltonian fromab initocalculations and performed tunneling calculations on it (using an early variation of instanton theory).While this work fell short of reproducing the experimental results (or resolving the discrepancy within them), it showed that there is a competition between the two mechanisms due to tunneling. Recently, thanks to the advance in computational power and path-integral based methods, on-the-flyab initioPIMD and ring-polymer instanton simulations have been applied to tackle this problem. PIMD shows that at high temperatures (i.e. 290 K and above), the stepwise pathway is favorable.[134,136]Over a wide range of temperatures, the two mechanisms form a tight competition.[136]As the temperature decreases,the contribution from the stepwise mechanism wanes, and at low temperatures (100 K and below), the concerted mechanism becomes dominate. The competition of the mechanisms could be why the DPT is difficult for experiments to characterize. It is worth pointing out that even at room temperature, both experiments and theory show that the DPT is enhanced significantly by quantum tunneling.[125,133,134,136]

The competition between concerted and stepwise transfer can be straightforwardly understood from the instantons.In panels a and b of Fig.3, minimum energy path(MEP)for the concertedtranstotranspath and the stepwisetranstocispath are shown,respectively,at 150 K and 100 K(black solid curves). The barrier for the concerted path is higher than that of the stepwise one. Stepwise path requires that the initialtransgeometry be thermally excited to an energy comparable to that of the intermediatecisgeometry for tunneling to take place. However, at low temperatures, concerted mechanism is dominant because there is not enough thermal energy available to reach the metastable IS via which the stepwise pathway goes through.

The main cause for qualitatively different DPT dynamics in prophyrin and porphycene originates from the distances between the nitrogens atoms in the inner ring, which are noticeably shorter in porphycene. The intramolecular H-bond strength are much stronger in porphycene than in porphyrin,resulting in qualitatively different PESs for DPT.Illustration of PESs(in a two-dimensional representation)for porphyrin and porphycene are shown in Fig.2. This indicates that it is possible to switch between stepwise and concerted mechanisms of proton/ hydrogen translocation by controlling the strength of hydrogen bonds between the reaction partners. This phenomenon could be exploited,e.g.,for the understanding of enzymatic reactions or in the design of proton wires.[131]

3.2. Porphycene on metal surfaces

STM is an important technique for studing the properties of the adsorbates. Molecules adsorbed on metal surfaces can be directly observed and monitored in STM experiments,therefore it is desirable to probe DPT in porphycene using STM.Due to the interaction between the molecule and the surface atoms, the properties of porphycene adsorbed on metal surfaces are different from that of gas- or condensed-phase porphycene.[48,50,138–142]For example,STM experiments and DFT calculations both show that thecisconfiguration is the most stable tautomer instead of thetransconfiguration on both Cu (110) and Ag (110) surfaces.[48,85,138]Also, porphycene adsorbed on the two surfaces have corrugated geometries that depends on the proton configuration instead of the planar geometry in the gas-phase.

Fig.4.Dependence of tautomerization rate on tunneling current,bias voltage,and temperature. (a) Trace of tip height (Z) measured for HH-porphycene.Traces are acquired at a set current (It) of 5 pA at Vbias =30 mV. (b)–(c)Schematic diagram of “high” and “l(fā)ow” states. (d) Trace of tip height (Z)measured for DD-porphycene. Right panels show magnified Z traces in the region marked by the dashed lines in the left,where a short-lived trans state is observed. Reprinted with permission from Ref.[48]. Copyright 2017 American Chemical Society,under the creative commons license 4.0.

Fig.5. Temperature (T) dependence of tautomerization rate for HH-, HD-,and DD-porphycene measured at It =10 pA and Vbias =10 mV. Reprinted with permission from Ref. [48]. Copyright 2017 American Chemical Society,under the creative commons license 4.0.

The rate of tautomerization of an adsorbed porphycene molecule can be directly monitored in STM experiments by tracing the change in the tunneling current or the change in the tip height over the molecule(since the molecule is slightly deformed from the planar geometry, as illustrated in Fig. 4(a)).For single porphycene molecule adsorbed on a Cu (110) surface,STM experiment showed that no spontaneous tautomerization(via tunneling)was observed at 5 K.Only when a relatively large bias voltage (190 mV) is applied, tautomerization starts to occur.[138]At higher temperatures(from 78 K to 95 K),thermally activated tautomerization can be observed in the STM experiment,and the observed rate show an Arrhenius behavior.[138]

In contrast to porphycene on Cu (110) surface, for porphycene adsorbed on Ag(110)surface,spontaneous tautomerization was observed directly via STM at low temperatures from 4 K to 20 K.The DPT process in this temperature range can be divided into two different regimes(classical-like regime and quantum-like regime): Above 10 K, the measured temperature dependence of the DPT rates shows an Arrheniuslike trend,suggesting that the underlying mechanism behaves classically. While below 10 K, DPT rates become almost temperature independent, indicating that the process is dominated by deep tunneling(Fig.5). In addition, a large KIE of about 100 between HH- and HD-porphycene have been observed, while the KIE between HD- and DD-porphycene is small (around a factor of 2).[48]These KIE factors corroborate to the Bell–Limbach model for the stepwise DPT mechanism. Furthermore, at 5 K, it has been directly observed that for DD-porphycene, the DPT process occurs via a stepwise mechanism. An additional level was observed in the tip height trace for DD-porphycene,which appears intermittently in between the“high”and the“l(fā)ow”states,as depicted in Fig.4(d).This state was attributed to the metastabletransconfiguration,which was not observed for HH-porphycene. This difference was attributed to more pronounced quantum tunneling in HHporphycene, which is believed to shorten the lifetime of the metastable intermediatetransstate.[48]

However, theoretical simulations reveal that DPT phenomena in porphycene adsorbed on Ag (110) are more complicated than previously known. Simulations with instanton theory clearly demonstrate a more competitive concerted tunneling mechanism under the lowTexperimental condition.[85]This indicates that similar to isolated porphycene molecule,there is also a competition between concerted and stepwise DPT at certain temperatures. Furthermore, instanton theory predicts that DPT mechanism can be divided into three different regimes at different temperatures, instead of just two regimes. As depicted in Fig.6,the numbers 1,2,and 3 illustrate the three different regimes as following:At high temperatures,stepwise transfer via classical hopping is dominant,displaying an Arrhenius behavior. As the temperature decreases,quantum tunneling becomes important, and the dominating mechanism changes to stepwise transfer via quantum tunneling. Despite that significant quantum tunneling occurs, the rate still display an Arrhenius behavior, albeit the slope differs from the previous regime. This is due to the metasable IS, as tunneling cannot break energy conservation, and thermal activation is still required to reach the IS. Further lowering the temperature, at cryogenic temperatures, concerted PT takes over as the dominate mechanism. This occurs because at cryogenic temperatures,there is not enough thermal energy available to reach the high-energy IS at which the stepwise path passes through.[85]However, discrepancies between experiment and theory remain. The calculated rates are much smaller than the measured rates in the experiments.[48]This discrepancy is very likely caused by the limitation of the DFT functional used in the calculations,as well as the effect of the STM tip on the system not being modeled.[85]Another discrepancy is that while the measured KIEs in the deep tunneling regime suggest that the underlying mechanism is stepwise (based on the Bell–Limbach model), the theoretically predicted DPT mechanism is concerted tunneling. We show later in this review that a similar discrepancy was found to exist in another system (PT in water tetramer adsorbed on an NaCl (001) surface). This indicates that either the Bell–Limbach model breaks down here or that there are other effects in the experiments missing in the theoretical works. For porphycene adsorbed on a Cu (110) surface, the calculated rates are in good agreement with experiment.[139]Both the calculated rates and the measured ones show an Arrhenius-like temperature dependence. Yet, theory revealed that under the experimental temperatures where the rates were measured, the underlying mechanism is stepwise tunneling,[85]rather than a classical hopping mechanism as usually indicated by the Arrhenius-like trend.

Fig.6. (a)Calculated instanton rates for the cis to cis DPT reaction of porphycene on Cu(110)between 75 K and 150 K.The inset shows experimental results from Ref.[138]. (b)Same as(a)but for Ag(110)between 7.5 K and 110 K.Orange circles represent the rate of the stepwise mechanism and blue squares of the concerted mechanism. The inset shows experimental results from Ref. [48]. (c) Schematic 1D potential energy surfaces of the stepwise(left) and concerted (right) reactions. The zero-point energy (ZPE) for reactant and products are shown. (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.Reprinted with permission from Ref.[85].Copyright 2020 American Physical Society,under the creative commons license 4.0.

3.3. DPT in tetraphenyl-porphyrin(TPP)on Ag(111)

TPP is a synthetic heterocyclic compound that resembles naturally occurring porphyrin.[143]Unlike porphyrin and porphycene, the energies of the twotransisomers of TPP differ, leading to an asymmetric barrier for DPT (Fig. 7). DPT process in TPP has been studied with STM for adsorbed TPP on Ag (111) (Fig. 8).[51]At low temperatures, the twotransisomers, referred to asHαandHκ, can be stably imaged and directly distinguished in high-resolution STM (Fig. 8). The conductance of the two isomers differ,therefore the tunneling current can be used to distinguish them,and the high and low conductance states correspond toHαandHκisomers respectively. When a bias voltage above a given threshold is applied,reversible switching between configurationsHαandHκcan be observed and traced by STM, indicating that current induced DPT occurs in this system.

Despite that thecisIS state was not observed in the STM experiment, the study claimed that DPT in TPP on Ag(111)occurs via a stepwise mechanism that goes through the metastablecisIS,which is similar to the DPT mechanism on porphyrin identified by previous experimental and theoretical studies.[114,115,144,145]They argued that judging from the previous experiments on prophyrin, thecisstate has a very short lifetime and relaxes either back to the initial state or to the final rotated configuration. Hence, the lifetime of thecisstate is too short to be observed in a typical STM experiment. It should be noted however, theoretical work for TPP adsorbed on metal surfaces is still absent, meaning that the mechanism of DPT has yet been fully verified. Future computational studies on TPP adsorbed on metal surfaces are desired to understand and provide more insights on the DPT mechanisms through asymmetric double-well potential.

Fig.7. The illustration of the stepwise DPT mechanism via cis IS in TPP(the phenyl groups are omitted for clarity). The macrocycle deformation on adsorption potentially lifts the degeneracy of both trans configurations.Reprinted with permission from Ref.[51]. Copyright 2012 Nature Publishing Group.

Fig.8. DPT in TPP on Ag (111). (a) STM image of configuration Hα. (b)Molecular model of configuration Hα. (c)STM image of configuration Hκ.(d)Molecular model of configuration Hκ. (e)Tunneling current versus time recorded at 21.9 V at the tip position indicated in panel (a). A switching between two current levels representing the high(h)and low(l)conductance states is clearly discernible.Reprinted with permission from Ref.[51].Copyright 2012 Nature Publishing Group.

4. Cyclic PT in systems without meta-stable IS

4.1. Water tetramer on Au supported bilayer NaCl film

The interaction between water and salt (NaCl) plays an crucial role in surface science,[146]environmental and water resources[147,148]and aqueous solution chemistry.[149]Thus, it is important to investigate PT phenomena in water tetramer adsorbed on NaCl surface. The concerted tunneling of four protons in water tetramer adsorbed on an Au-supported NaCl(001)film has been directly visualized by state-of-the-art STM technique.[49]Two chirality states of the water tetramer exist,i.e.,the clockwise state(CS)and the anti-clockwise state(AS),as depicted in Fig.9. According to experiment,the Clterminated tip is positioned slightly off the center of the water tetramer. The switching dynamics of water tetramer chirality was monitored by recording the tunneling current as a function of time. When the tip is far above the water tetramer,the current remains constant,indicating no switching between the two chiral states. Once the tip height is reduced, two current levels appear, in addition to a sudden increase in the tunneling current. The low and high current levels can be assigned to the CS and the AS, respectively (as shown in Fig. 9(b)).Then,The switching rates between them can be measured,detailed information please see Ref. [49]. The measured rates are weaklyT-dependent below 20 K,which was attributed to deep quantum tunneling in the PT process.

To understand whether the chirality switching of tetramer happen through concerted PT, stepwise PT, or rotation of the water molecules, DFT calculations have been carrier out on the system.[49]The potential energy profile of these switching mechanisms has been calculated, with the STM tip also modeled in the calculations. It is found that the barrier height(～0.8 eV)is its lowest when the four protons concertedly hop from the H-bond donor to the acceptor. While for stepwise PT, the static energy profile is depicted in Fig. 10. Starting from the stable non-defective tetramer, PT along one hydrogen bond would result in a paired defect state(one hydroxide defect and one adjacent hydronium defect). However,this reaction path has a significantly larger barrier height (～2 eV)compared with that of the concerted PT path,which suggests a very small possibility of tetramer switching via stepwise PT either classically or quantum mechanically. Besides concerted and stepwise PT,rotation of four H-bonds has also been studied. It has been noticed that although the barrier height of rotation pathway is comparable to that of the concerted PT pathway, however, the former has a larger effective barrier width,which is not ideal for tunneling to occur(details can be found in Ref.[49]). Also,the water rotational barrier is considerably more sensitive to the tip height compared to the PT process,whereas the experiments found that the rate is rather insensitive to it. Therefore,it was concluded that concerted PT is the most likely mechanism for the chirality switching in the water tetramer observed in the experiments.

Fig.9. (a) Schematic showing manipulation of the chirality of the tetramer by a Cl-terminated tip. Left: the tetramer stays in the clockwise state (CS)when the tip is far away from the tetramer.Middle:reducing the tip height by 230 pm leads to chirality switching. Right: lifting the tip back to the initial height leaves the tetramer in the anticlockwise state(AS).(b)Tunneling current trace recorded during the chirality manipulation shown in panel(a). Two current levels can be clearly distinguished, where the low and high current levels correspond to CS and AS,respectively. Left and right insets: Adsorption configuration(upper)and STM images(lower)of CS and AS tetramers,respectively.The green stars in the STM images denote the tip position where the current trace is acquired. O, H, Au, Cl-, and Na+ are denoted by red,white, golden, cyan, and blue spheres, respectively. Reprinted with permission from Ref.[49]. Copyright 2015 Nature Publishing Group.

Fig.10. Potential energy profile along the stepwise PT pathway from CS and AS state.In a stepwise PT process there are two paired defect states(PD1 and PD2)and a separated paired state(SPD).Insets show the structures of these states. Reprinted with permission from Ref. [49]. Copyright 2015 Nature Publishing Group.

Yet, further studies on the KIEs of this system raises questions concerning the underlying mechanism. It has been found that the chirality switching rate of the water tetramer reduced by over an order of magnitude after replacing one H2O molecule with D2O.Further when all four H2O molecules are replaced by D2O,the switching rate does not further decrease(Fig.11).This observed behavior resembles the Bell–Limbach model for stepwise PT, which is add odds with the previous conclusion that the underlying mechanism is concerted PT. This unexpected observation has drawn new attention to this system. Theoretical calculations on the rates and KIEs of the chirality switching in the water tetramer on NaCl system has been performed usingab initoPIMD and instanton theory.[150,151]Quantum mechanical switching rates obtained from PIMD calculations(QTST rates to be precise[96])at different temperatures down to 50 K, show that strong quantum tunneling occurs in this process, even at room temperatures.KIEs are computed at 50 K, showing that although the rate ratio for 4H2O and 3H2O+1D2O,k4H/k3H1Dis close to the experimental result,[150]the rate for 4D2O(k4D)is significantly lower thank3H1D. KIEs predicted by instanton theory agrees with this finding as well,[151]and the theoretical KIEs qualitatively agree with the Bell–Limbach model for concerted PT(Eq.(5)).

Fig.11. Comparison of the quantum and classical chirality switching rates for water tetramer on NaCl over a range of temperatures down to 50 K.The label “PIMD” refers to QTST rates,[96] which are RPMD rates without the recrossing factor, computed from quantum free energies obtained from constrained PIMD simulations. They can be viewed as upper limits for the RPMD rates of this process. Red and black dashed lines represent the calculated classical and quantum switching rates of 4H2O tetramer, respectively. Blue up-triangle and green down-triangle are calculated switching rates for 3H2O+D2O and 4D2O tetramers. Reprinted with permission from Ref.[150]. Copyright 2018 AIP Publishing.

Some possible sources of the discrepancy between experiment and theory have been discussed. Since thermal rate theories could not reach agreement with experimental findings, it might happen that microcanonical rates are actually measured in the experiments instead of thermal rates. This is a possibility due to the extremely low temperature in the experimental set up and excitations from the electric current applied in the system for probing. Microcanonical instanton theory[152,153]has been used to compute the rate as a function of excitation energy instead of temperature. We can obviously see from Fig. 12(b) that at any excitation energy, the rate for 3H2O+D2O is close to the rate for 4H2O,while being much higher than the rate for 4D2O. This trend is qualitative the same as the trend in thermal rates calculated by PIMD and instanton theory and also in agreement with the Bell–Limbach model. Therefore,the discrepancy between theory and experiments unlikely result from energy excitation caused by the STM tip.[151]Previous studies also considered other important factors, for example, the influence of the height of STM tip,electric field between the tip and the sample,underestimation of the reaction barrier of generalized gradient approximation (GGA) exchange-correlation interactions, etc, and concluded that these factors are unlikely to have any qualitative impact.[49,150]

Fig.12. (a)Canonical instanton rates. The temperature dependence of chirality switching rates. Due to numerical difficulties, the point marked with an open symbol could not reach our convergence criteria and is shown only for the purpose of a qualitative reference. (b)Microcanonical instanton rates.The red,blue,and green solid lines show the rates of chirality switching for 4H2O,3H2O+D2O,and 4D2O,respectively.Reprinted with permission from Ref.[151]. Copyright 2022 AIP Publishing.

Meanwhile, some speculations on the causes of the discrepancy were also made.[151]For example, since the experiment is at very low temperatures(5 K),if the coupling between the water tetramer and the substrate phonons is weak enough,it is possible that the tunneling phenomenon is coherent.[128]In this case, the observed rate could display different behaviors and could be influenced by the measurement technique.One thing certain is that more efforts from both the theoretical and experimental sides are desired in the future.

4.2. PT in hexagonal ice

Hexagonal ice(iceIh)is the most common phase of water at low temperature and ambient pressure,constituting a nearly perfect tetrahedral network of hydrogen bonds. The protons in iceIhare disordered,there are six possible orientations for a water molecule in its tetrahedral bonding environment,each corresponding to a different arrangement of protons in the four H bonds.[154]IceIhprovides a rich platform for the study of NQEs due to its peculiar arrangement of H-bonds governed by the ice rules.

An incoherent quasielastic neutron scattering (QENS)measurement on iceIhreported the existence of unexpected non-harmonic motion of hydrogen at low temperatures in iceIh.[7]A significant broadening of QENS signal in an iceIhsample was observed at 5 K, while such broadening is absent in the QENS signal of a partially deuterated sample at the same temperature(Fig.13). They proposed that this anomaly is caused by concerted PT of 6 protons in the hexagonal ring,as illustrated in the inset of Fig.13,and that partial deuteration breaks the symmetry as well as the concerted PT. They also estimated that this process has a large rate of 2.7×1011s-1,and affects only a few percent of the total number of hydrogen atoms.[7]This result was astonishing, as prior studies never found significant proton mobility in ice at low temperatures and ambient pressure.

Fig.13. Quasielastic contribution in hydrogenated ice Ih (°)compared with partially deuterated(◇)at 5 K.Inset: a sketch of the proposed concerted tunneling of the hydrogen atoms in the ordered loops.Reprinted with permission from Ref.[7]. Copyright 2009 American Physical Society.

The first-time reporting of an experimental observation what appears to be rapid concerted PT in ice soon attracted enormous attention in the field. A theoretical study invested this system usingab initioMD and PIMD simulations.[19]They computed the classical and quantum free energy barriers for the PT process in iceIh, along a PT collective variable defined using the average position of the 6 protons in the H-bonds (Fig. 14). At 50 K, the quantum free energy barrier reduces significantly to only a quarter of the classical free energy barrier. The ring-polymer stretches across the barrier in the simulation, indicating the dominance of deep tunneling. Even at the room temperature and above, NQEs still play a role (albeit limited role), as the quantum free energy barrier is slightly reduced compared to the classical one.Importantly, a comparison between PIMD simulations in the low-temperature tunneling regime and in the high-temperature classical regime reveals qualitative differences in the transition mechanism. At high temperatures, the preferred mechanism is “stepwise”, with decoupled single-proton hopping sequentially along individual H-bonds,creating charge defects on the hexamer ring. However, at low temperatures, the 6 protons tunnel concertedly, largely circumventing charge defects and the associated energy penalty,making it the favorable mechanism at these temperatures.

Fig.14. Free energies as a function of the global collective variable Φc from quantum(blue squares: 50 K;red triangles: 320 K)and classical(black circles: 300 K, green diamonds: 50 K) simulations, detailed definition of Φc could be found in Ref. [7] Fig. 1; linear dotted lines are guides to the eye.The left and right insets schematically show the proton configuration of the ordered six ring in the. Reprinted with permission from Ref.[19]. Copyright 2014 American Physical Society.

PIMD simulations have also been used to investigate H/D isotopic-substitution effects on the mechanism of the concerted tunneling of six protons in proton-disordered iceIh.[6,155]Free energy profile of fully protonated iceIh(H system) partially deuterated (H/D system) has been computed.Results show that partial deuteration of the H system leads to inequivalent protons,thereby suppressing concerted proton tunneling. If one H on the water hexamer ring is substituted with D,charge separations starts to appear in the hexamer during the PT process(Fig.15(a)),leading to a larger free energy barrier compared to the case with 6H(Fig.15(d)). These findings,as well as a few other theoretical works,[156,157]support the experimental observation[7]that concerted PT occurs in iceIh.

Fig.15. (a)–(c) Representative transition-state configurations of the 5H1D(a),6D(b),and 6H(c)systems. Oxygen and hydrogen atoms[(a)–(c)],and formal H3O+ and OH- ionic defects [(a) and (b)] are shown in red, gray,blue, and green, respectively. The overlaid spheres represent quantum delocalization according to the usual path-integral representation[(a)and(c)].(d)Free energy profiles along the global proton-transfer coordinate Δ at 50 K(see Ref. [155] for the definition), for 5H1D: red triangle; 6D: black circle;6H: blue diamond. Adapted with permission from Ref. [155]. Copyright 2014 John Wiley and Sons.

Althoughab-initioPIMD simulations clearly show that concerted transfer becomes favorable in iceIhat low temperature,the enormously large rate claimed in the experiment has yet been explained. The PIMD calculations for concerted PT in ice[6,19,155]computed only free energy profiles,but did not provide any estimate on the rate of the concerted PT process.Therefore, it still remains an open question on whether the concerted PT happens in ice can occur on a time scale observable in experiments. Especially since the early QENS experiment claimed an extraordinarily high rate for the concerted PT process,[7]the controversy around this topic has not ended.A recent QENS experiment by Kolesnikovet al. revisited this system,which found that the earlier QENS experiment by Boveet al.[7]cannot be reproduced,and that there is no concerted proton tunneling in iceIhat cryogenic temperatures.[57]State-of-the-art QENS measurements(with the energy resolution similar to and four times better than the energy resolution in the original experiment of Bove and co-workers)show that at 1.8 K and 5 K, QENS broadening cannot be observed in iceIh. In addition, the measured spectra for the iceIhand the reference vanadium sample were almost identical,indicating that no proton tunneling in iceIhat temperatures down to 1.8 K. Besides QENS broadening and the measured spectra,the heat-capacity is another measurable quantity that can reveal information on tunneling.Possible presence of the tunneling splittings in iceIhshould be reflected in the heat-capacity dataCP(T)at low temperatures. If the concerted tunneling of six protons exist through a symmetric well,the tunneling splitting of the ground state,described as a two-level energy state(TLE) with energy splitting ΔE, would be discernible from anomalies in the heat capacityCP(T)at low temperature. For a TLE with energy splitting ΔE, the tunneling splitting contributes to heat capacity at low temperature as

whereNTLE= 0.04 is a number of TLE states andkBis a Boltzmann constant. The measured low temperatureCP(T)(from 0.5 K to 10 K) for iceIh[158]shows no anomaly at all.It is clear thatCP(T)follows the dependenceT3,as it should be for dielectric material with Debye behavior[159]for vibrational density of states at low energies. Figure 16 shows comparison of the experimentalCP(T)for iceIh[158]with the calculated contribution of TLE states if the tunneling of water protons existed. ΔEvalues of 0.4, 0.2, and 0.062 meV were used. It is apparent that the experimental data certainly disagree with the calculatedCTIE(T). Thus, one can concluded that there is no proton tunneling based on the heat capacity at low temperatures,[57]adding another piece of evidence against the claim of concerted PT in iceIh.

Fig.16. Heat capacity of ice Ih at low temperatures(blue points and curve)are from Ref.[158]and calculated contribution to CP(T)due to presence of 4% of two-level states (Schottky anomaly) with splitting energy. Reprinted with permission from Ref.[57]. Copyright 2018 American Physical Society.

5. Summary and outlook

This review covered several representative example systems for(cyclic)multiple PT,with the major focus on the theoretical modeling of these processes,meanwhile covering key experimental findings. A major point of interest in (cyclic)multiple PT processes is the understanding and characterization of the underlying mechanism. There are two types of PT mechanisms in these systems, stepwise PT and concerted PT,and possible competition between the two,which can give rise to rich behaviors and phenomena. Theoretical insights on these systems as well as accurate modeling of multiple PT mechanisms are not only crucial to the development of theoretical tools with high predictability, but also important for the designing of nanostructures with functionality(i.e.molecular switches,proton conduction wires). Thanks to the developments of both experimental techniques and theories (especially path-integral based methods), in the past decade or so,some significant advances in the characterization and understanding of cyclic PT mechanisms have been made. For gasphase porphyrin and porphycene,experiments have found that the PT mechanism in the two differs qualitatively despite the two molecules being similar,with the former showing a clear stepwise PT while the latter is difficult to firmly characterize.PT in the two systems can be modeled accurately with theory,showing that indeed the stepwise mechanism is favorable in porphyrin,while the two mechanism form a close competition on porphycene due to quantum tunneling. Porphycene on metal surfaces can also be modeled quite well with theory(e.g.instanton theory), with the predicted PT rates on Cu (110) in good agreement with STM experiments. On Ag(110),theory predicted an interesting change of mechanism with temperature,showing three different regions.

Proton switching in water tetramer on NaCl is another interesting system that attracted great attention in recent years.While initially both experiments and theory show that the process precedes via a concerted mechanism, later, KIEs were probed experimentally,and the discrepancy between the measured and predicted KIEs raises new questions on the mechanism. To date, this is still an open question (mystery) surrounding this system. Another question that is still under debate is on whether concerted PT occurs in a hexagonal ring in iceIh. An early experiment claimed observations(indirectly)of such events, and that theory also found that such process is indeed possible. However, later experiments could not be reproduced. From the theory perspective, while there is no doubt that concerted PT can occur in ice, the rate which has not been evaluated,leaving questions regarding the timescale of such event. These open questions regarding multiple PT mechanism requires more efforts in the future from both the experimental and theoretical side to resolve. It is worth noting that apart from the development of better quantum rate theories in the future, the development of more accurate and efficient electronic structure methods (including improvements to DFT functionals) that can be applied to metallic and condensed phase systems is also crucial.

Acknowledgments

Project supported by the National Basic Research Programs of China (Grant No. 2021YFA1400503), the National Natural Science Foundation of China (Grant No. 11934003),the Beijing Natural Science Foundation(Grant No.Z200004),and the Strategic Priority Research Program of the Chinese Academy of Sciences(Grant No.XDB33010400).