Anomalous rotational-symmetry breaking in proton arrangement of surface-confined cyclic hydrogen bonds revealed by atomic force spectroscopy

Domain: physics
Published: 2025-12-02
Editorial outcome: major_revision

Abstract
In hydrogen-bonded materials and biosystems, microscopic proton ordering, which is strongly influenced by the quantum nature of protons, underpins diverse macroscopic phenomena, including phase transitions, chemical reactions, biomolecular processes and collective proton transfer. Yet resolving proton arrangements and characterizing their quantum behavior at the atomic scale remain challenging due to the small size of protons and the lack of effective approaches. Here, we exploit bond-resolved atomic force microscopy and spectroscopy (BR-AFM/AFS) to probe signatures of proton ordering in surface-confined benzimidazole (BI) assemblies. By performing BR-AFS along the apparent H-bond between proton donor and acceptor nitrogen atoms, we extract information consistent with H-bonding directionality and signatures compatible with quantum proton delocalization. We observe anomalous rotational-symmetry breaking in the proton order of the cyclic hexamers, arising from the coexistence of both localized and quantum-delocalized protons. The chirality of a single hexamer can be reversibly switched by altering its adsorption registry combined with collective transfer of six protons. Path-integral molecular dynamics calculations unravel that nuclear quantum effects promote proton delocalization and concomitant reduction of donor-acceptor distance. Our findings enable atomic-level detection of complex proton orders, engineering proton-based quantum states, and elucidation of long-range proton transfer mechanisms in molecular architectures.

Introduction
Hydrogen bonds are profoundly influenced by nuclear quantum effects (NQEs) owing to the low mass of the proton1–3. To fully understand many fundamental biological and chemical processes mediated by H-bonds–such as enzymatic reactions4–6, signaling7, and the dynamics of water8–the quantum nature of protons must be considered. In two- or three-dimensional molecular crystals, where intermolecular H-bonding motifs create specific proton arrangements9–12, NQEs can significantly impact proton ordering, thereby dramatically altering structural and physical properties. For example, quantum delocalization drives proton order-disorder transition and H-bond symmetrization in high-pressure phases of ice13, and promotes high-temperature superconductivity in sulphur hydrides14. Recent high-resolution scanning probe studies of small molecules on atomically flat surfaces in ultrahigh vacuum (UHV) have provided preliminary evidences for proton sharing, including water-hydroxyl complexes on Cu(110)15, long-range-ordered Zundel cations on Au(111) and Pt(111)16, and H-bonded one-dimensional 2,5-diamino-1,4-benzoquinonediimine chains on Au(111)17. However, these observations have so far relied primarily on imaging, while complementary spectroscopic evidence remains limited.

Another challenge is to extend the investigation of proton ordering and NQEs from small interfacial molecular systems, such as water18, to those with heavier skeletons10. In this regard, the imidazole ring represents an ideal candidate. Capable of both donating and accepting a proton, it serves as an important building block in biological systems, facilitating low-barrier H-bonds (LBHB)4,5. Moreover, single-component crystals of imidazole derivatives can achieve above-room-temperature ferroelectricity, with electric polarization switchable via cooperative proton transfer through the H-bonds (N−H···N⇄N···H−N)19.

In this work, we introduce an approach based on bond-resolved atomic force microscopy and spectroscopy (BR-AFM/AFS) for determining H-bond directionality in benzimidazole (BI) architectures on Ag(111). We reveal the conventional unidirectional H-bonding in linear BI chains, while in cyclic hexamers, we unravel ground-state proton orders featuring coexisting localized and quantum-delocalized protons (states 1 and 4 in Fig. 1). We further demonstrate that a single hexamer can interconvert reversibly among four proton-ordering states (Fig. 1). Changing the adsorption registry suppresses proton delocalization, whereas a collective six-proton transfer not only reverses the chirality of the cyclic H-bonds but also restores the non-classical H-bonds (Fig. 1). Our work highlights the important role of NQEs and the local environment on proton ordering, as well as correlated nature of many-body proton transfer20,21 within interfacial molecular architectures.

Results
Probing H-bond directionality via force spectroscopy
The five-membered imidazole ring in the BI molecule contains two nitrogen atoms. For clarity, we refer to the pyridine-like nitrogen as “N site” and the pyrrole-like nitrogen as “NH site” (see Fig. 2a). Upon adsorption on Ag(111), the prochiral BI molecule transforms into a chiral entity, denoted as either the l- or r-enantiomer (Fig. 2a). At submonolayer coverage, BI molecules self-assemble into various structures, including extended molecular chains, closed loops and six-membered clusters (Fig. 2b), driven by attractive intermolecular N−H···N hydrogen bonding provided by the imidazole ring22. The chains grow along the 112 direction of Ag(111) (Fig. 2b), indicating a preferred epitaxy. To specify the N−H···N bond directionality in this system, it is essential to distinguish between the N and NH sites in the imidazole ring. Fig. 2c depicts a constant-current STM image of a single BI molecule recorded with a CO-functionalized tungsten tip (see also Supplementary Fig. 3), permitting submolecular resolution23. Note that a single-lobe feature appears on both sides of the molecule, adjacent to the N and the NH sites, with the latter being more prominent (cf. the molecular structure in Fig. 2c and the detailed discussion below). Using frequency shift (Δf) imaging in constant-height mode24, both pentagonal and hexagonal rings are clearly resolved (Fig. 2d), with slight expansion and distortion due to CO bending at the tip apex25,26. However, the N and NH sites cannot be differentiated.

To distinguish between pyridine and pyrrole-like nitrogen atoms, we performed site-specific force-distance spectroscopy [FS or F(z)], which records Δf changes as a function of tip-sample distance27–30. Fig. 2e depicts a set of F(z) spectra recorded along a path that traverses two corner sites of the imidazole ring (see dashed arrow in Fig. 2d). The left projection in Fig. 2e displays the Fmin values extracted from each F(z) spectrum as a function of the lateral tip position d along the path. The Fmin(d) curve exhibits an asymmetric M-shape, with two peaks of unequal height separated by ≈4.6 Å, matching the apparent distance between two nitrogen sites in the imidazole ring (cf. Fig. 2d). The right projection in Fig. 2e shows two F(z) curves extracted at the maxima of the two peaks in the Fmin(d) curve (cf. blue and red points in left projection in Fig. 2e), with Fmin(N) − Fmin(NH) ≈ 7 pN (cf. Supplementary Fig. 4). FS provides contrast consistent with differences between N and NH sites mainly through atomic-scale electrostatic interactions with the CO tip, reflecting the distinct electron density distributions of these two sites, as elaborated below. Notably, our results are consistent with prior FS characterization of naphthalocyanine molecule, which distinguishes the inner nitrogen sites with and without hydrogen31.

Next, we applied line-FS to determine the H-bonding directionality in an open chain (Fig. 2f and Supplementary Fig. 5). Although the bright intermolecular contrast in BR-AFM image cannot be directly interpreted as an H-bond32–34, the probe-particle (pp) model reveals that the “bright line” arises from the relaxation of the probe particle (CO molecule) above the potential landscape formed by two neighboring atoms26. When the sharp line coincides with an H-bond, it traces the path between the proton donor and acceptor atoms. Accordingly, FS measurements were performed along the bright connections between neighboring imidazole rings (Fig. 2f, g), crossing two adjacent nitrogen atomic sites. Along the apparent H-bonds, characteristic M-shaped Fmin(d) curves are also obtained with the first peak consistently higher than the second (left projection in Fig. 2h; Supplementary Fig. 6d). Two F(z) curves recorded at two peak positions show a clear difference at Fmin (right projection in Fig. 2h), with Fmin(N) − Fmin(NH) ≈ 4 pN. Therefore, the H-bond directionality in the BI chain adopts a left-pointing N···H−N configuration (Fig. 2g and Supplementary Fig. 6). This assignment is corroborated by two additional observations: (i) the appearance of a dark depression at the right-side termination, attributed to a pyridine-like N atom35 (white arrow in Fig. 2f); (ii) the apparent length of the single (C−N) bond being greater than that of the double (C=N) bond (Supplementary Fig. 7)25. Although the experimentally measured apparent distance between intra- or intermolecular N and NH sites are systematically larger or smaller than the DFT-calculated values (Fig. 2a and Supplementary Fig. 6e) due to CO-tip bending, this does not hinder the accurate discrimination between the two sites.

Anomalous symmetry breaking in BI-hexamer proton order
Next, we focus on self-assembled BI hexamers, which are homochiral owing to their H-bonding motif. These hexamers adopt four different adsorption configurations (Supplementary Figs. 8 and 9), denoted as ∓8.5°-L/R or ∓40°-L/R. The rotation angle (∓8.5° or ∓40°) is defined between the high-symmetry axes of the hexamer and the close-packed directions of the Ag(111) surface (Fig. 3b; Supplementary Fig. 8), while the clockwise (+) or counterclockwise (–) turn is associated with handedness (R/L) of the cyclic H-bonds in the hexamer36,37. Fig. 3a, b show a BR-AFM image and the DFT-optimized adsorption geometry of the −8.5°-L hexamer. The center of the hexamer occupies an hcp hollow site (cf. Supplementary Fig. 10), giving rise to threefold rotational symmetry: the benzene ring in L1, L3 and L5 sits over a hollow site, whereas in L2, L4 and L6 it sits on a top site. Notably, the +8.5°-R also adopts an hcp-centered geometry with C3 symmetry (Supplementary Fig. 10b), and DFT calculations show that both enantiomorphic configurations are planar on the surface with nearly identical adsorption energies (Supplementary Fig. 11).

To determine the cyclic H-bonding directionality in the hexamer, we first analyzed bond-order contrast in experimentally obtained BR-AFM images but found no systematic trend (Supplementary Figs. 12 and 13). By contrast, applying the same analysis to AFM images simulated with different effective tip charge38 reveals a clear bond-order contrast consistent with the expected chiral motif (Supplementary Fig. 14), and we clarify this discrepancy below. To gain further insight, we performed line-FS characterizations (see also Supplementary Figs. 15–18). Fig. 3g depicts six datasets, measured consecutively for the six H-bonds in the −8.5°-L following a clockwise sequence (red arrow in Fig. 3a). Each panel displays a characteristic M-shaped Fmin(d) curve along with two F(z) spectra recorded at two peak positions (red and blue points and curves in Fig. 3g). An overall trend can be discerned: along the measurement path, the second peak in the M-shaped Fmin(d) curves is predominantly higher than the first. Therefore, the Fmin(d) curves show features compatible with a clockwise N−H···N bonding scheme in the −8.5°-L (Fig. 3b). For comparison, line-FS measurements were also conducted for the +8.5°-R hexamer, revealing a counterclockwise N−H···N bonding scheme (Supplementary Fig. 15d, h).

To further support our interpretation, we calculated electrostatic potential (ESP) of the −8.5°-L hexamer on Ag(111) using DFT and simulated both AFM images and site-dependent FS with the PPAFM code39 and a negatively charged CO tip40,41. The ESP map reveals that the N site has a more negative potential than the NH site (Fig. 3e). In particular, the simulated M-shaped Δfmin(d) curve with unequal peak heights for all six H-bonds (Fig. 3h, i), together with the corresponding Δf(z) curves at the N and NH sites (Fig. 3j), reproduces the experimental results (Fig. 3g), thus validating the line-FS approach for determining H-bond directionality.

A closer inspection of Fig. 3g shows that, for H-bonds L3-L4 and L4-L5 in the −8.5°-L hexamer, the two Fmin(d) peaks have nearly equal heights: |Fmin(N) − Fmin(NH)| ≲ 1 pN, (Supplementary Figs. 15i and 17c; Supplementary Note 2). For such H-bonds, the extracted F(z) curves at peak positions overlap near their minima (orange-highlighted plots in Fig. 3g). At a tip-sample distance corresponding to Fmin, the AFM image of −8.5˚-L shows a single-lobe feature corresponding to the N-site in the imidazole ring on the right-hand side of the BI molecule for H-bonds L1-L2, L5-L6 and L6-L1. By contrast, a double-lobe feature is clearly observed in L3-L4 and L4-L5 and faintly in L2-L3 (Fig. 3c). Similar patterns are also observed for the +8.5˚-R hexamer (Supplementary Fig. 15f). However, the equal-peak M-shaped Fmin(d) curves and the corresponding double-lobe AFM features deviate from the DFT-based PPAFM simulations (Fig. 3f, h). Such observations imply that the donor and acceptor nitrogen sites become indistinguishable, reminiscent of both symmetric H-bonds15,42 and LBHBs4,5, and suggest the possibility of proton delocalization within the H-bond. AFM imaging near Fmin, together with line-FS characterization on multiple as-grown ∓8.5°-L/R and ∓40°-L/R hexamers, consistently reveals anomalous rotational-symmetry breaking of proton order (Supplementary Figs. 19-21). This contrasts with the nearly C6-symmetric, flower-like pattern from the PPAFM simulation (Fig. 3f) and with the computed M-shaped Δfmin(d) profiles with unequal peak heights expected for equivalent H-bonds (Fig. 3h). Moreover, neither small biases (Supplementary Fig. 22) nor a CO-terminated tip41 perturb the broken-symmetry proton order, although it is not discernible by tunneling spectroscopy (Supplementary Fig. 23). In different hexamers, one, two or three of the six H-bonds exhibit signatures of proton delocalization. Across a total of 56 cases, the mean fraction of non-classical H-bonds is ≈40% (see Supplementary Note 3).

Unravelling proton delocalization via PIMD simulations
To examine the influence of NQEs on proton ordering in the BI hexamer, we performed ab initio molecular dynamics simulations using both conventional and path-integral schemes (AIMD/PIMD)43 at 50 K. Compared to AIMD, which considers nuclei as classical particles, PIMD treats both electrons and nuclei quantum mechanically, providing a more accurate description of proton behaviours17,44. To simplify the system, we modelled a gas-phase BI dimer (L3-L4) with its initial geometry extracted from the optimized −8.5°-L-hexamer/Ag(111) (Supplementary Fig. 24a; see also “Methods”). Fig. 4a depicts the probability distributions of a proton between two nitrogen atoms–N (donor) and N′ (acceptor)–as a function of reaction coordinates RC derived from PIMD and AIMD calculations. The RC is defined as RC = (RNH − RN′H) / (RNH + RN′H) (1) where RNH and RN′H are the instantaneous NH and N′H distances, respectively (see inset in Fig. 4a). The PIMD calculations reveal that the proton is delocalized between two nitrogen atoms, in stark contrast to the AIMD results, where the proton remains covalently bonded to the donor atom (Fig. 4a, b). Free energy is calculated from ΔF = −kBT ln P(RC) (2) where P(RC) is the probability distribution of RC and kB is the Boltzmann constant. Under full quantum description, the free energy profile features a double-well valley at RC ≈ −0.18 and 0.16, separated by a small free energy barrier of ≈ 6 meV at RC = 0 (Fig. 4b). This may considerably enhance proton tunneling, leading to the occurrence of proton delocalization (see a snapshot from the PIMD simulations in Fig. 4a). Conversely, the AIMD-calculated free-energy profile displays a single valley at RC ≈ −0.18, indicating that thermally activated proton transfer is rare. This is consistent with the relatively large classical reaction energy barrier (≈ 32 meV) computed for proton transfer in the H-bond of the BI dimer using the climbing image nudged elastic band (CI-NEB) method with zero-point energy (ZPE) correction (Supplementary Fig. 24b).

We further analyzed the temporal variations of NN′ distance and RC (Fig. 4c) to reveal the connection between NN′-distance distribution and proton delocalization. At the quantum mechanical level with PIMD, the averaged NN′ distance decreases moderately from ≈ 2.7 to ≈ 2.6 Å, accompanied by the occurrence of proton delocalization, where the averaged RC approaches zero (orange-highlighted region in Fig. 4c). For AIMD calculations, the NN′ distance oscillates around an equilibrium value of ≈ 2.7 Å (upper plot in Supplementary Fig. 24c). Intriguingly, the time-dependent RC evolution obtained from quantum simulations reveals a change in the sign (lower plot in Fig. 4c), corresponding to a proton tunneling process from N to N′, which is absent in the AIMD calculations (lower plot in Supplementary Fig. 24c). Within the timescale (10−13 s) of the quantum simulation, proton tunneling from N to N′ only occurs once due to the low free energy barrier (≈ 6 meV), which corresponds to a transfer rate on the scale of 1012 s−1 based on the quantum transition state theory method applicable to water tetramer45. Accordingly, we infer a rapid back-and-forth dynamic behavior of delocalized protons, which underlies the nearly overlapped F(z) spectra recorded at two nitrogen sites of the non-classical H-bonds (orange-highlighted plots in Fig. 3g and Supplementary Fig. 15). Since line-FS integrates over time, it reflects the time-averaged proton positions.

Notably, in BR-AFM images non-classical H-bonds appear shorter than classical ones (Supplementary Fig. 25), consistent with PIMD findings (Fig. 4c) and prior theoretical predictions44. However, such geometrical effects cannot be reproduced at the DFT level (Supplementary Fig. 26). To capture the coexistence of localized and delocalized protons in cyclic H-bonds (Fig. 3c, g), we performed PIMD simulations of a constrained BI-hexamer/Ag(111) surface at 5 K using the reaction force field46 (Methods; Supplementary Fig. 27). Starting from the DFT-relaxed −8.5°-L hexamer (Fig. 3b), the N···N distances of L3–L4 and L4–L5 were shortened slightly relative to other N···N spacings (Fig. 4d; and Supplementary Fig. 27a). The resulting proton probability distributions and free-energy profiles for L3–L4 and L4–L5 are presented in Fig. 4e, f. First, proton delocalization occurs only in shortened H-bonds (≈ 2.4 Å)42, whereas for N···N ≳ 2.8 Å the proton remains localized at the donor atom (Fig. 4d). Second, L3–L4 and L4–L5 display unequal proton probability distributions and corresponding asymmetric free-energy profiles. In the gas-phase dimer (Fig. 4a), removing one proton yields two mirror-symmetric configurations that interconvert via proton tunneling, producing a symmetric probability distribution. In the hexamer, however, L5 retains a localized proton on its left side, whereas L3 is proton-deficient on its right side. Proton tunneling could therefore transiently generate a benzimidazolium cation at L5 and a radical BI species at L3. Although such fleeting states are not resolved in time-averaged FS measurements, they may account for the suppressed proton delocalization in L3–L4 and L4–L5.

Controllable transitions among four proton-ordering states
We demonstrate that a single BI hexamer can reversibly interconvert among four distinct proton-ordering states: two ground states containing non-classical H-bonds and two intermediate states in which classical H-bonds prevail (Figs. 1, 5). At liquid nitrogen temperature, STM scans induce conformational change of hexamers between two equivalent adsorption configurations (Supplementary Fig. 8, and Supplementary Movie 1). This implies that a single hexamer can switch between two enantiomorphic structures when the chirality of the cyclic H-bonds is reversed through collective transfer of all six protons (cf. Supplementary Fig. 11). To verify this assumption, we performed tip manipulations on hexamers at 5 K (Supplementary Figs. 28, 29; and Supplementary Movie 2). Fig. 5a presents an AFM image of a −8.5˚-L hexamer in its ground state (state 1; cf. Fig. 3c). After tip manipulation, the hexamer rotated 17° around its center (Fig. 5b). Line-FS analysis of this configuration (+8.5°-L) reveals that the clockwise H-bonding directionality persists, while non-classical H-bonds are significantly suppressed in this intermediate state 2, as evidenced by the downshift of the ΔFmin values for L3-L4 and L4-L5 relative to state 1 (shaded area in Fig. 5h). Intriguingly, +8.5°-L can be transformed into +8.5°-R (state 4) via multiple proton transfer induced by inelastic tunneling electrons47. The CO-tip was positioned above the NH site of a BI molecule (green dot in Fig. 5b), and the bias voltage was ramped from zero to 0.5 V while keeping the feedback-loop open. A sudden jump in the Δf(V) curves was observed at 0.36 V (Fig. 5g and Supplementary Fig. 30), matching the NH stretching mode energy (Supplementary Fig. 31). This operation re-established the non-classical H-bonds, yielding +8.5°-R ground state (Fig. 5c and Supplementary Fig. 32b, f). The reverse transition (state 4 to state 1) follows a similar pathway (Figs. 1, 5d–f; and Supplementary Figs. 33, 34), proceeding via intermediate state 3. Importantly, multiple-proton transfer could not be triggered in states 1 and 4, indicating that ground-state proton ordering of hexamer is interlocked with its adsorption configuration.

To rationalize why classical H-bonds dominate in states 2 and 3 (Fig. 5h, i), we performed additional DFT calculations (Supplementary Fig. 35). We found that, in states 2 and 3, three pyridine-like nitrogen atoms in the hexamer occupy top sites on Ag(111) (highlighted in yellow in Supplementary Figs. 32c, 33c), whereas in states 1 and 4, all pyridine-like nitrogen atoms reside at the hollow sites (Supplementary Fig. 11). This indicates a stronger interaction between the proton-acceptor nitrogen atom and the substrate in the intermediate state configuration, which is corroborated by electron-density-difference analysis, Δρ = ρAg+hexamer − ρAg − ρhexamer (3) where ρ represents electron density of each component. Δρ cross sections along the H-bond (∓8.5°-L) indicate greater electron-density redistribution beneath the top-site proton-acceptor nitrogen atom (N′) than that beneath its hollow-site counterpart (cf. dashed rectangles in Supplementary Fig. 35b, d). Strong electronic hybridization causes a pronounced tilting of the N′ atom’s σ orbital (lone pair) out of the in-plane geometry (a black arrow in Supplementary Fig. 35d), thus disfavoring proton tunneling.

Methods
Sample preparation
Ag(111)/mica substrates were employed and the Ag(111) surface was prepared by repeated Ar+ ion sputtering at 0.9 keV and annealing at ≈ 700 K for several cycles. Benzimidazole (Sigma Aldrich, 98%) was fully degassed and then evaporated at around 320 K onto the substrate maintained at room temperature. At low coverage, the molecules favored the formation of hexamers, whereas high coverage led to the formation of extended chains and closed rings.

Synthesis of 1H-benzo[d]imidazole-1-d
The synthesis of 1H-benzo[d]imidazole-1-d (deuterated at the N1 position) involves exchanging the N–H proton with deuterium using a suitable deuterium source. To achieve selective N1-deuteration of unsubstituted benzimidazole, we employed an acid/base-controlled, divergent deuteration protocol specific to N-unsubstituted imidazoles, using trifluoroacetic acid-d (TFA-d) as the deuterium source (Supplementary Figs. 1, 2)51. This method enables selective labeling at the N1 position without affecting the aromatic rings or the C2 position. The reaction was carried out under acidic conditions with CF₃CO₂D (60 equivalents) at 130°C for 12 h. The extent of deuteration was confirmed by ¹H NMR spectroscopy, where the disappearance of the N–H proton signal indicated successful incorporation of deuterium (for details see Supplementary Note 1).

STM/AFM measurements
The scanning tunneling microscopy and atomic force microscopy (STM/AFM) measurements were conducted in a Polar STM/AFM system (Scienta Omicron) at 4.8 K with a base pressure of better than 2 × 10−10 mbar. The qPlus sensor used had a resonance frequency (f0) of 24.5 kHz, a spring constant (k0) of ≈ 1800 Nm−1, and a quality factor (Q) greater than 150000 at 4.8 K. A netched tungsten (W) tip was used in the experiment, which was pulsed and poked into the metal surface to ensure a clean, symmetric, and atomically sharp front edge. Oscillation amplitudes of 0.5 Å were employed. All STM and AFM images were recorded in constant-current and constant-height modes, respectively, using a CO-terminated tip unless otherwise stated. The bias voltage was applied to the sample with respect to the tip. AFM images were acquired at a bias voltage of 0 V. The scanning height (z) of AFM images is relative to STM setpoint on the bare Ag(111) surface, with z > 0 indicating an increase in the tip-sample distance, and z < 0 indicating a decrease in the tip-sample distance. Image processing was carried out by using Nanotec WSxM software52.

Statistics and reproducibility
All STM and AFM measurements were repeated on multiple samples and with different tips, yielding consistent results. The representative images shown in the main text were selected from reproducible datasets obtained under similar experimental conditions.

BR-AFS measurements
Force spectroscopies [FS or F(z)] were measured with Us = 0 V. To perform line-FS on an H-bond, a BR-AFM image was first obtained for the H-bond. Following this, FS measurements were conducted along the bright line between the two molecules. Finally, another bond-resolved AFM image was recorded after line-FS to ensure no noticeable drift had occurred. Short-range F(z) were obtained by subtracting long-range FS recorded on Ag(111) and converting Δf(z) to F(z) via the Sader-Jarvis method53.

Rotation and collective proton transfer via tip manipulation
Before manipulation, the tip was cleaned/sharpened by repeated pulsing/poking on Ag(111). For an as-grown hexamer (∓8.5°-L/R), the tip was placed over one molecule at Us = 10 mV, It = 10 pA; feedback was then disabled and the tip lowered until It ≈ 20 nA to enlarge tip-molecule interaction. Lateral motion at 1 nm s−1 rotated the hexamer into an intermediate state (±8.5°-L/R), which preserves cyclic H-bond chirality. We then imaged the hexamer, positioned the tip at an NH site, and ramped the bias until an abrupt Δf change at |Us| ≈ 0.35 V indicating collective proton transfer; a follow-up AFM image confirmed the ground state (±8.5°-R/L). The protocol was repeated > 30 times with high reproducibility.

DFT calculations and AIMD simulations
DFT calculations were performed using the Vienna ab initio simulation package (VASP)54,55, with the projector augmented wave (PAW) pseudopotentials56,57. The exchange-correlation effect was described by the Perdew-Burke-Ernzerhof exchange-correlation functional58 within the generalized gradient approximation (GGA). The modelled adsorption structure consists of a BI hexamer adsorbed on the three-layer Ag(111) substrate of a 7√3 × 12 unit cell and the thickness of the vacuum slab is larger than 11 Å. During structural optimization, the substrate was fixed, and the cutoff energy of the plane wave basis set was set to 520 eV. Due to the large scale of the absorption model, the structural optimizations were considered converged when the maximum force on each atom was less than 0.01 eV Å−1. Only the Γ point in the Brillouin zone was sampled. DFT-D3 method with Becke-Johnson damping function59 was used to account for the van der Waals interaction for the surface adsorption.

In ab-initio molecular dynamics (AIMD) based on the classical nuclear approximation and the climbing image nudged elastic band (CI-NEB) calculations, the cutoff energy was set to 600 eV and the degrees of freedom in z-direction were fixed. The maximum force on each atom was less than 0.005 eV Å−1 after the CI-NEB calculation had reached convergence. The AIMD simulations were performed for 1935 fs with a timestep of 0.1935 fs at a temperature of 50 K. The canonical (NVT) ensemble was sampled using the Langevin thermostat60,61.

Ab-initio PIMD calculation
Ab-initio based PIMD simulations employing 8 beads were performed using the i-PI software62, combined with the time-dependent ab initio package (TDAP)63–65 implemented in Quantum Espresso66,67. The SG15 Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials68,69 was used to describe the core electrons. The Perdew-Burke-Ernzerhof functional58 was adopted to describe the exchange-correlation effect. The kinetic energy cutoff for wavefunctions and for charge densities were set to 100 Ry and 400 Ry, respectively. The PIMD simulations were performed for 385 fs with a timestep of 0.1935 fs at 50 K. The canonical (NVT) ensemble was sampled with the Path Integral Langevin Equation (PILE) thermostat70. PIMD based on the force field was performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)71 and the reaction force field46.

The BI dimer was derived from the optimized −8.5°-L structure and pre-optimized in the gas phase prior to AIMD/PIMD simulations. In the calculations, one proton was removed from the BI dimer, resulting in a negatively charged species. This modification places the remaining proton in an asymmetric environment, which promotes quantum delocalization within the intermolecular hydrogen bond. At each timestep, the positions of the N, N′, and H atoms were extracted, the RC = (dNH − dN′H) / (dNH + dN′H) was determined, and all possible combinations were enumerated. The resulting probability distribution and free energy profiles of the reaction coordinate were smoothed using a 3-point discrete convolution kernel with a weight of [0.2, 0.6, 0.2] across 200 intervals, followed by a sliding rectangular window filter spanning 9 bins.

PIMD simulation based on force field
Force field-based PIMD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)71,72 and the reaction force field46. The normal mode path-integral driver73 implemented within LAMMPS were used to carry out the simulations. A Nose-Hoover massive chain thermostat74 was used to sample the canonical (NVT) ensemble of 32 beads at a temperature of 5 K. The total duration of the simulations was 20 ps with a timestep of 0.2 fs. During the simulations, the molecular backbone (C and N atoms) and the substrate (Ag atoms) were fixed. The probability distribution and free energy profiles of the reaction coordinate were smoothed by applying a 3-point discrete convolution kernel with a weight of [0.1, 0.8, 0.1] across 200 intervals. The simulated BI-hexamer on Ag(111) surface exhibited an asymmetric structure (Fig. 4d and Supplementary Fig. 27), where the threefold rotational symmetry was broken by rotating molecules L2, L3, L5, and L6 by ≈ 4.8° towards L4 around the center of the hexamer.

AFM simulations
AFM images and force–distance spectra were simulated with the open-source probe-particle model26,39. Parameters: the effective lateral stiffness k = 0.50 Nm−1, radial stiffness = 30 Nm−1; oscillation amplitude A = 100 pm. The CO-functionalized tip was modeled as a dz2-like charge quadrupole41 with Q = −0.10 e (where e denotes the elementary charge). Tip–sample interactions combined Lennard–Jones potentials with DFT-based electrostatics computed for BI hexamer on Ag(111).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The paper investigates hydrogen bonds in cyclic hexamers of benzimidazole (BI) molecules on a Ag(111) surface using\nSTM and AFM with CO functionalized tips. Resolving the directionality and discerning classical hydrogen bonds (N−H···N)\nfrom non-classical hydrogen bonds with quantum proton delocalization (N···H···N) is claimed. Manipulation between different\nhydrogen bonding motifs by atom manipulation is investigated.\nThe quality of the data is very good. However, the interpretation of the results and the claimes made are not very convincing.\nI do not think that there is enough experimental evidence for the distinction between classical and non-classical hydrogen\nbonds. It would be great, but is not convincingly shown. Several questions remain open and need to be considered. I\nrecommend that the paper is reconsidered after mayor revision.\nOverall, the data is of high quality. However, I find the assignment of classical hydrogen bonds (N−H···N) and non-classical\nhydrogen bonds with quantum proton delocalization (N···H···N) not convincing. I also see a potential problem with hexamers,\nthat show both type of hydrogen bonds. If molecules have a classical hydrogen bond on one side, and thus either have a\nhydrogen or no hydrogen bonded to that N atom, it requires the other N atom of that molecule to have no or a hydrogen,\nrespectively, and thus also a classical hydrogen bond, or they would become radicals, right? I can imagine a structure with\nall non-classical hydrogen bonds. However, in the mixed hexamer, with both types of bonds, I see a conceptual problem.\nThis needs to be discussed in more detail, or maybe reconsidered.\n1. I think it is also very instructive to consider the closed shell Kekule structure of benzimidazole. In it, the two bonds that\nconnect the carbon atom at the 5 membered ring to N and NH site are different. A formal double bond to the N and a formal\nsingle bond to the NH group. The resulting bond-order related contrast with AFM is in line with the AFM image of the chain\nshown in Fig. 2f. The formal double bonds appear shorter compared to the single bonds. Moreover, the related difference in\nelectron density above the two N atoms might contribute to the contrast mechanism on N and NH sites, as observed in the\nmonomer. The bond order of N-C bonds in the hexamers could be analysed.\n2. In the hexamer the differences in bond-order contrast are less compared to the open chain. This should be discussed.\nMight it imply that possibly in all hexamer structures all the H bonds are delocalized?\n3. In a non-classical hydrogen bond with quantum proton delocalization (N···H···N), would it need to be symmetric? Could\nnot the proton be delocalized, but still be (in average) closer to one of the two N, allowing some asymmetry of the N···H···N\nbond?\n4. A structure that has in part classical hydrogen bonds (N−H···N) and in part non-classical hydrogen bonds with quantum\nproton delocalization (N···H···N) as claimed in states 1 and 4 in Fig. 1 seems problematic to me. The molecules connecting\nnon-classical and classical H bonded segments have a “wrong” count in hydrogen atoms. That is having one classical\nhydrogen bond in the ring should determine all other classical hydrogen bonds in the ring. It seems to me that non-classical\nhydrogen bonds with quantum proton delocalization (N···H···N) should only be possible for either all bonds in the hexamer,\nor none. Please discuss in detail how a molecule that on one side has a classical H bond, can have a delocalize one on its\nother side.\n5. Related to the previous point: In theory the hydrogen bonds with quantum proton delocalization (N···H···N) is only obtained\nif one hydrogen is removed from one molecule (Fig. 4b). It seems this “trick” is needed because of the problem described in\nthe previous point. What is the charge state of the structure in Fig. 4b? Is it a neutral doublet, or charged?\n6. Could a hydrogen possibly indeed be removed from a hexamer? Providing another possible explanation for the observed\nswitching and asymmetric contrast in the hexamers (and relate to the different 8.5 and 40 deg tilted structures).\n7. The authors make not clear how they distinguish between classical hydrogen bonds (N−H···N) and non-classical\nhydrogen bonds with quantum proton delocalization (N···H···N). More specific: How do the authors come to the specific value\nof 1 pN, that they use to discern classical and non-classical N bonds? Why this specific value? Possibly, (N···H···N) bonds\nmight be asymmetric. Moreover, the molecules might be tilted slightly on the surface or their local adsorption environment\nmight be different possibly contributing to the Fmin difference, that is not related to the nature of the hydrogen bond. The\nvalue of 1 pN to discern seems chosen arbitrary.\n8. The authors state: “Notably, we found that the centre of the +8.5°-R also coincides with a hcp site (Fig. 3d and\nSupplementary Fig. 7). This indicates that the mirror symmetry for freestanding 6-membered enantiomorphs does not hold\nwhen considering their adsorption registries (Fig. 3c,d).” I do not understand the argument. If the center would not be the\nsame site, I would agree that the structures can not be mirror symmetric. I do not see why both structures being centred\naround hcp sites makes the authors conclude they are not mirror symmetric. In fact they do look mirror symmetric to me. (Flip\nat horizontal line in Fig. S7).\n9. Related to the previous point: Is the energy difference calculated for the enantiomers in Fig. S8 significant? Or is it within\nthe noise (energy cutoff criterium) of the calculation?\n10. The energies for adsorption under +/-8.5 deg (Fig. S8) could be compared with similar calculations for the local\nadsorption minima expected at +/-40 deg.\n11. The authors write: “Thus, our methodology demonstrates that H-bond directionality can be unambiguously determined\nvia line-FS measurements, offering a clear resolution to the longstanding debate over whether the bright intermolecular\ncontrast represents an H-bond or an artifact due to CO-tip flexibility [31,36].” First of all: what is the resolution to the debate?\n(I fear the authors might mean that the contrast is due to the hydrogen bond). What is the reasoning? In this context, please\nconsider that you also state: “This feature was observed in the Fmin(d) curves measured for a single BI molecule and linear\nchain (Fig. 2e,h), where Fmin(NH) is lower than Fmin(N).” Note, that you observed the feature also in a BI monomer, that is\nin the absence of any hydrogen bonds.\n12. Is the adsorption site of the structures shown and switched in Fig. 5 identical? That is, are the transition states 2 and 3\nalso centered around hcp sites?\n13. Related to previous point: Can the 40 deg oriented structures also be switched or manipulated to become other\nstructures, e.g. the +/-8.5 deg structures?\n14. If orbital resonances could be accessed by STM, might that provide additional experimental data on the type of hydrogen\nbond?\n15. Optional: The paper could very likely be greatly improved by performing experiments with deuterated molecules. Effects\nof proton tunneling could be revealed very convincingly by the comparison of proton vs. deuteron tunneling. It could help\ndiscerning conventional from delocalized hydrogen bonds. In the present manuscript the assignment and distinction of these\ntwo bonds is not clear.", "reviewer_2": "The authors present an experimental and theoretical investigation of quantum proton ordering in hydrogen bonds using\nbond-resolved atomic force microscopy and spectroscopy (BR-AFM/AFS). They demonstrate that enantiomorphic\nbenzimidazole (BI) hexamers exhibit unique quantum proton ordering in their ground states, which can be manipulated via\nscanning probe techniques. These experimental findings are supported by ab initio molecular dynamics (AIMD) and path-\nintegral molecular dynamics (PIMD) simulations, confirming the presence of nuclear quantum effects (NQEs) in the system.\nThis work is highly relevant to the study of quantum effects in hydrogen bonds and molecular self-assembly on surfaces. The\nmanuscript provides an important insight into quantum proton ordering and hydrogen bonding at the atomic scale. However,\nto ensure the validity of the conclusions, the authors should address major concerns. Additionally, minor corrections and\nrefinements in figures and terminology are necessary.\nMajor Comments\n1. The primary issue with the manuscript is the interpretation of hydrogen bond directionality based on force spectroscopy.\nAs the authors are likely aware, there is ongoing debate in the community regarding whether the contrast observed above\nhydrogen bonds originates from the bonds themselves or other effects. Specifically, it remains unclear whether CO-tip AFM\ntruly images hydrogen bonds or if the observed contrast stems from other intermolecular interactions\n(see, e.g., S. K. Hämäläinen et al., Phys. Rev. Lett. 113, 186102 (2014)). Since this assumption is fundamental to the manuscript's\nconclusions, the lack of supporting data beyond point force spectroscopy is surprising. A straightforward way to validate\nthese findings would be to provide a series of BR-AFM images at different tip-sample distances to corroborate the\nspectroscopy results. Additionally, the Probe Particle Model is a well-established method for simulating BR-AFM images\nwith high accuracy. Given that this code is open-source, the authors could perform their own simulations to support their\ninterpretations. Interestingly, my own simulations of BI in the gas phase did not reveal significant differences between the N\nand N-H sites, suggesting that the observed contrast might be a topographical effect rather than a direct measure of\nhydrogen bonding. Further scrutiny of the provided images suggests that the molecule may not be lying flat on the surface,\nwhich could influence the total force measurements. Another noteworthy observation is the discrepancy in the minimum\nforce values extracted from Figures 2e and 2h, which range from ~10 pN to ~40 pN. The authors should clarify the physical\norigin of this variation—could it be due to incomplete background signal removal? Additionally, why should the N-H site be\nmore attractive to the CO tip than the N site? This could be better understood by analyzing the electrostatic potential (ESP)\nof the molecule or performing density functional theory (DFT) calculations of the force spectroscopy. The authors discuss the\nrole of adsorption registry in proton ordering, yet the specific electronic interactions between the substrate and the molecules\nremain unclear. Does the Ag(111) surface actively stabilize particular proton configurations? Are charge transfer effects\ninvolved in influencing hydrogen bonding? Additional discussion or computational analysis would strengthen these claims.\n2. The manuscript suggests that proton transfer in BI hexamers is primarily induced by inelastic electron excitation from the\nSTM tip. However, alternative mechanisms, such as thermally activated or substrate-mediated processes, cannot be ruled\nout. Have the authors conducted control experiments at varying temperatures or tip biases to exclude these possibilities?.\nCan the authors also provide a mechanism explaining how inelastic electron tunneling induces\nproton transfer?\n3. Since all experiments are conducted at 5 K, it is unclear how temperature affects proton delocalization and tunneling.\nWould quantum proton ordering persist at higher temperatures, or would classical hydrogen bonding become dominant?\nA\ndiscussion on this aspect is necessary.\nMinor Comments\n• Figures 3 and 4 require improved annotations to clearly indicate key measurements, particularly the extracted force-\ndistance parameters.\n• The phrase “broken rotational symmetry” appears frequently but is not well defined. The authors should clarify its precise\nmeaning and implications.\n• The introduction should provide a more detailed comparison with previous studies on quantum proton delocalization in\nsimilar molecular systems.\n• Typographical and Formatting Issues:\no Page 3, line 7: “proton delocalisation” → “proton delocalization” (for consistency with American English).\no Page 7, line 23: “image” → “images”; “Large” should be lowercase.\no Page 7, lines 24-25: The word “products” appears redundantly.\no Figure 3F appears misaligned compared to Figure 3B—this should be corrected.\no Figure 5: The red arrow is not explained in the caption.\no Supporting Information: Some materials (e.g., Scheme 1, Figures S1 and S2) are not referenced\nin the main text. Their relevance should be clarified.", "reviewer_3": "The study presents a compelling investigation into quantum proton ordering in hydrogen-bonded benzimidazole hexamers\nusing bond-resolved atomic force microscopy and atomic force spectroscopy. The combination of high-resolution\nexperimental techniques and path-integral molecular dynamics (PIMD) simulations provides strong evidence for proton\ndelocalization and nuclear quantum effects in a well-defined surface-confined system. A major strength lies in the ability to\nmanipulate and switch between quantum and classical proton-ordering states using an AFM tip, offering insights into proton\ntransport and potential applications in quantum materials. In my view, this research makes a significant contribution to the\nunderstanding of molecular-scale quantum effects and merits publication, as it introduces both novel insights and technical\nadvancements in the field. However, I have several questions and comments that should be addressed before acceptance,\nwhich I list below.\nThe study is restricted to benzimidazole hexamers on Ag(111), which may not fully represent biological or bulk hydrogen-\nbonded systems. Surface interactions with Ag(111) may influence proton behavior in ways that do not generalize to solution-\nbased or gas-phase environments. I would appreciate a comment from authors whether the observed quantum proton\nordering extends beyond the Ag(111) surface? Would the same proton ordering be observed on weakly vs. strongly\ninteracting surfaces? Also, would this behavior persist in a freestanding benzimidazole hexamer or in bulk crystalline\nenvironments? In any case, I believe a short justification on the choice of Ag(111) surface would be appreciated by the\nreader.\nThe interaction of the CO-functionalized tip with the system could alter proton behavior, raising questions about the true\nground-state configuration in the absence of an external probe. Could authors comment on that?\nHow does the temperature as well as external electric field affect proton ordering? Could an external field stabilize or disrupt\nthe quantum proton sites? In case of temperature I would appreciate a comment why temperature of 50K was chosen for\nsimulations.\nHave authors consider using other techniques to further confirm that the observation is truly quantum mechanical behavior?\nFor example, inelastic neutron scattering or high-resolution vibrational spectroscopy could confirm the quantum vs. classical\nnature of the discussed proton behavior.\nI also wonder whether an isotope effect would influence the findings? Would for example deuteration lead to observable\nshifts in proton ordering, verifying the role of quantum effects? Could authors comment on that?\nOn a smaller note. Perhaps I am exaggerating, but in Fig. 4b and 4c, I notice that the PIMD simulation results appear to show\ntwo maxima/minima on the acceptor side. I am curious whether this is simply an artifact of the simulation or if it could indicate\nadditional effects that have not been discussed in the main text. Could this possibly be related to proton delocalization\nasymmetry or subtle interactions within the H-bond network? I also strongly suggest to mark in these figures as well as Fig.\n4d which side corresponds to the donor and which to the acceptor, as this would help in correctly interpreting the data\nwithout looking into them main text.\nIn Fig. 3e,f, I find it somewhat challenging to clearly distinguish the single and double lobes that the authors refer to in the\ndiscussion. While Fig. 3i,j certainly help clarify the interpretation by quantifying the force differences, perhaps there is a more\neffective way to visually enhance these features in the AFM images.\nThe following sentence is missing a verb: “It is important to note that although the experimentally determined apparent\ndistance between intra- or intermolecular N and NH sites systematically larger or smaller than the DFT-calculated values…”\nIs there a reason the authors use \"Supplementary Fig. 3\" instead of a shorter format like \"S3\"?\nIn the Supplementary Information, the authors use hcp in italics, but in the main text, the formatting is inconsistent. It would\nbe best to standardize the notation throughout the paper.\nIt may just be my impression, but in Fig. 2e, it looks like the N site and NH site are not exactly chosen at the maxima, but\nrather slightly next to them. Is this just an optical illusion, or could the authors clarify how these positions are precisely\ndefined? Additionally, I noticed a difference in the donor-acceptor distance when comparing Fig. 2e and Fig. S2c. in Fig. 2e,\nthe value is 460 pm, while in Fig. S2c, the value is 490 pm. Is this difference correct, and if so, what is the reason behind it?", "reviewer_4": "The manuscript titled \"Detection and control of quantum proton ordering in hydrogen bonds at the atomic scale\" presents a\nhighly complete and data-solid study that utilizes bond-resolved atomic force microscopy and spectroscopy (BR-AFM/AFS)\nto explore quantum proton ordering in hydrogen-bonded systems. In my opinion, the experimental data are of high quality,\nand the observed phenomena are both novel and fascinating. The theoretical calculations, particularly the path-integral\nmolecular dynamics (PIMD) simulations, provide robust support for the experimental findings, offering a detailed mechanistic\nunderstanding of proton delocalization and quantum effects in hydrogen bonds. Overall, this is an excellent piece of work,\nand I agree with its publication in Nature Communications. However, several issues need to be addressed to optimize the\nmanuscript further:\n1.Regarding the long-standing debate over the origin of the clear resolution in CO-tip-based chemical bond detection (i.e.,\nwhether the bright intermolecular contrast represents a hydrogen bond or an artifact caused by CO-tip flexibility), the\nmanuscript mentions that \"line spectroscopy measurements can unambiguously determine hydrogen bond directionality,\nproviding a clear resolution to this debate.\" Could the author comment more on this conclusion? Why can local force\nspectroscopy be more convincing in determining the position of the H atom than the features in constant height imaging like\nFig. 2f?\n2.his work reports a remarkable observation of symmetry-broken quantum proton ordering in the ground states of chiral\ncyclic H-bonds within BI hexamers on Ag(111), despite the symmetric adsorption registry as shown in Fig. 3c-d. To clarify\nthis intriguing result, could the authors comment more on this phenomenon? Could it be a local distortion of this hexamer\nduring the deposition of the molecules or a local imperfection of the Ag surface?\n3.The PIMD simulations in Fig. 4 elegantly demonstrate proton delocalization in a gas-phase BI dimer (L3–L4). However,\nthe quantum behavior of the full BI hexamer adsorbed on Ag(111)—where six H-bonds interact cooperatively—might not be\nfully captured by this simplified dimer model. Could the authors comment more on whether the dimer-based simulations are\naccurate enough to describe the cooperative transfer of protons in a hexamer?\n4.The experimental observations of symmetry-broken H-bonds and their manipulation (e.g., Fig. 3, 5) were conducted at 4.8\nK, while the PIMD simulations (Fig. 4) were performed at 50 K. Why did the author choose this simulation temperature and\nwhether the choice of simulation temperature (50 K) impact the agreement between the predicted proton delocalization\nbehavior and the experimental results obtained at ultralow temperatures (4.8 K)?\n5.The manuscript highlights the critical role of adsorption geometry in stabilizing non-classical H-bonds (e.g., symmetric H-\nbonds in hexamers vs. classical H-bonds in chains; Supplementary Figs. 12–13). Could tailored modifications of adsorption\nsites (e.g., substrate step edges, strain engineering) enable the emergence of symmetric H-bonds even in chain\nconfigurations? And why is the delocalized proton only found in the hexamer?"}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors provided extensive new data and explanations and clarified most points.\nI think the addition of the bond-order related contrast of C=N and C-N bonds improved the MS.\nThe DFT calculations showing a flat-lying geometry of the hexamers (Supplementary Figs. 10 and 17), is an important\nimprovement.\nThe most important point in my opinion is the assignment of the switching to proton tunneling. It is great that the authors\nperformed experiments with deuterated molecules. Showing a different voltage threshold for H and D tunneling would be\nconvincing proof of the switching being related to proton/deuterium tunneling. This the authors attempt to show in the new\nFig. S34. However, with the few individual reported voltages for switching and the way they are presented it is not\nconvincing. This could be improved. It is not clear if that is a proof or not. I think without taking additional data this question\ncan and should be answered before publication. It would be nice and a very convincing proof of proton/deuterium tunneling\nif different voltages thresholds could indeed be extracted from the experiment. But, I would also not object to publication, if\nnot.\nMy points 1-4 refer all to Fig. S34:\n1. The voltages of switching should better (or at least in addition) be shown against the current at which the switch occurred.\nAssuming a process mediated by tunneling electrons, the relation to the tunneling current would be more direct than to the\ntip heigh. It would also remove possible effects that might be related to the set point being on a different surface (Ag) and at a\ndifferent voltage compared to the voltage at which the reaction occured with the tip on the molecule, and most important\neffects due to different tips being used. Please plot switching Voltages against the current at which the switch occurred.\nAs a side note: The different used tip heights z, z’, Z make it extremely hard to understand and add confusion. The tip height,\nz’ seems given wrt. z with seems given wrt. the STM setpoint – one has to find different passages of the MS and SI to find that\nout. Is this here z, or z’ or Z? And what are the differences?\n2. Assuming an electron mediated process (which the authors do, and I agree with that), likely with a small yield, the voltage\nof switching will be stochastic. That is, the individual switching voltage will vary in different individual experiments (therefore\nI would not name that individual switching voltage “voltage threshold”), even under identical experimental conditions.\nMoreover, the switching voltage will even depend on the ramp velocity of the voltage. That means:\n2a) The (statistical error of the points in Suppl. Fig. 34e) should be estimated. Ideally this is done by repeating switching with\nidentical parameters several times (at least for one condition). I assume that authors have such data. Please provide error\nbars for the points in Supplementary Fig. 34e (and plot against current).\n2b) The parameters of the bias sweeps (including the time for the sweep) should be similar or at best identical for the\nexperiments on H and D molecules. And the parameters should be provided.\nI assume the trend in Figure 34e will be the same as now when plotting threshold Voltages as a function of current. Which\nwould be more convincing. However, not sure how large the error bars will be. This cannot be inferred with the information I\nhave. I assume the error bars will be of the same order of magnitude as the observed differences in switching voltages\nbetween H and D molecules.\nMinor points in Suppl. Fig. 34.\n3. Should not the labels “NH” better read “ND”, as these experiments were on deuterated molecules?\n4. On the y-axis in panel e, the values for 280 and 320 mV seem swapped.\nOther comments:\n5. It seems at least three different tip heights, z, z’ and Z in Suppl. Fig. 15 in some other instances (Suppl. Fig. 34) are used.\nThe use of the different z, and z’ (and Z) seems extremely inconsistent and unclear to me. I understood from the text z’ is\nused for spectroscopy, but it is then mentioned as tip height for many AFM images in the SI. E.g. in S. Fig. 9, z is used in the\npanels, but z’ in the caption. Both seem to denote the same thing and have the same values. But in the in the method section\nit says that z and z’ differ by 30 pm. Please check all instances and make clear what is the difference between z and z’.\nDo you really need both of them, why not subtract the 30 pm and only use z, in all instances?\n6. I focused on Fig. S34 and reporting the things to be improved only in tis Figure took me a while. I did not do that so\ncarefully on all figures. Please check yourself. One other thing I did spot: Fig. S7: Delta f scales in c and d are very different\ncompared to the other panels. Please double check, if correct. If correct, maybe comment on it. Was a different tip used?\n7. I now realized that the method to distinguish N and NH sites is very similar to what was shown in F. Mohn, L. Gross, N.\nMoll, and G. Meyer, “Imaging the charge distribution within a single molecule,” Nat. Nano. 7, 227–231, 2012, see\nSupplementary Figure S2. Compare e.g. with Figure S15e. This seems to support the analysis of the current MS. The work\nshould be cited and mentioned in this context. I think it also strengthens the assignment. Note that in the cited work, the\npositions of N and NH groups could be inferred independently by orbital density maps measured with STM. Moreover, the\nrelation of the AFM contrast, to the KPFM contrast (mainly related to electrostatic interaction) was measured and discussed.", "reviewer_2": "The authors have substantially revised the manuscript, adding height-dependent AFM/FS data, isotope substitution,\nextended simulations (DFT + PPAFM, PIMD, classical MD), and additional controls. These additions significantly improve\nthe robustness of the study.\nThe central claim—that AFM line frequency-shift profiles can reveal proton delocalization driven by nuclear quantum effects\n—remains controversial, as electrostatics, adsorption registry, and fast hopping are not entirely excluded. Nevertheless, the\nexpanded dataset now provides a consistent and technically rigorous case that certain AFM signatures are consistent with\nnon-classical hydrogen bonds, and that these signatures cannot be captured by classical models or simple geometric\nvariations. With careful framing, this contribution will be of broad interdisciplinary interest and is suitable for publication in\nNature Communications.\nI therefore recommend acceptance after minor revision, contingent on the following clarifications and additions:\nGeneral:\nThe manuscript establishes a correlation between AFM/FS signatures and the presence of non-classical H-bonds. However,\nstating that AFM directly resolves or discerns proton delocalization seems too strong. I encourage the authors to moderate\nthe wording, using terms such as ‘consistent with’ or ‘compatible with’ rather than categorical claims.\nSections where the language should be softened\nAbstract (p.2): “…BR-AFM/AFS to directly probe the proton ordering…”\nSuggest: “…BR-AFM/AFS to probe signatures of proton ordering…”\nAbstract (p.2): “…we determine H-bonding directionality and detect signature of quantum proton delocalization.”\nSuggest: “…we extract information consistent with H-bonding directionality and signatures compatible with quantum proton\ndelocalization.”\nResults, “Elucidating hydrogen bond directionality…” (p.4): “…FS distinguishes N from NH sites…”\nSuggest: “…FS provides contrast consistent with differences between N and NH sites…”\nResults, hexamer analysis (p.6):\n“…Fmin(d) curves reveal a clockwise N−H···N bond directionality…”\nSuggest: “…Fmin(d) curves show contrast compatible with a clockwise N−H···N bonding scheme…”\nSummary (p.10):\n“…a line-FS based method that resolves H-bonding directionality…”\nSuggest: “…a line-FS based method that provides evidence consistent with H-bonding directionality…”\nSome concersn that still need clarification from the authors:\nQuantitative robustness of the classification criterion:\nReport the number of spectra averaged per bond and provide error bars for |ΔF_min|.\nSummarize the fraction of bonds classified as non-classical across all measured.\nHeight-dependent analysis:\nPresent Δf(z) traces at N vs NH with residuals from your fits, to show explicitly the z-range over which the classification is\nstable. Indicate the uncertainty from the Sader–Jarvis conversion.\nInternal charge effects:\nInternal charge redistribution within a single BI molecule is expected to distort AFM contrast patterns directly impacting into\nbond-length analysis from experimental images. Did the authors explicitly check for such distortions and compare with bond-\nlength analysis obtained from calculations? Please clarify.", "reviewer_3": "The authors have satisfactorily addressed the points I raised in my previous report. As I am not a specialist in the specific\nexperimental methodology, I am not in a position to comment in detail on some of the technical issues discussed in\nresponse to other reviewers’ comments and will leave their evaluation to experts in the field. From my side, I have no further\ncomments.", "reviewer_4": "The authors have thoroughly addressed all comments raised in this round of review. They provided further/newer data\nwhere necessary and added relevant explanations where required. The paper reads more clearly, and the analysis of the\nresults is more convincing. This study advances our understanding of the properties of cyclic hydrogen bonds. I therefore\nrecommend publication of “Anomalous symmetry breaking in proton ordering of surface-confined cyclic hydrogen bonds” in\nNature Communications."}, "author_response": "Point-by-point response to the reviewers’ comments (NCOMMS-24-85724-T)\nWe sincerely thank all the reviewers for their valuable time and positive evaluation of our work.\nTheir insightful comments and suggestions have helped us improve the clarity and scientific rigor\nof the manuscript significantly. We address all raised questions by providing additional\nexperimental and computational data, along with further analysis. Figure numbering refer to those\nin the revised manuscript and supplementary information (SI).\nNote: In the first section titled “Responses to Common Questions”, we address several shared\nquestions that have been slightly rephrased and reorganized for clarity.\nIndex:\nResponse to Common Questions: (pages 2-5)\nPoint-by-point response to the Reviewer #1’s comments (pages 6-13)\nPoint-by-point response to the Reviewer #2’s comments (pages 13-17)\nPoint-by-point response to the Reviewer #3’s comments (pages 17-22)\nPoint-by-point response to the Reviewer #4’s comments (pages 22-24)\nCommon question 1: Why can hydrogen bond directionality be determined by force spectroscopy?\n(Reviewers #1, #2 and #4)\nReply: The determination of hydrogen bond (H-bond) directionality relies on distinguishing\nbetween the N and NH site in the imidazole ring. Our reasoning is as follows:\n(1) We agree that the bright intermolecular contrast observed in CO-tip AFM images cannot be\ndirectly ascribed to hydrogen bonds (H-bonds). We are aware that similar bright-line features\nalso appear in systems without H-bonds, such as silicon chains on Ag(110)1.\n(2) Importantly, the probe-particle (pp) model reveals that the bright intermolecular contrast arises\nfrom the relaxation of the probe particle (CO molecule) above the potential energy landscape\nformed by two neighboring atoms2. When the sharp line’s position coincides with an H-bond,\nit traces the path between the proton donor and acceptor atoms.\n(3) Pyridine-like (N) and pyrrole-like (NH) nitrogen atoms can be distinguished using force\nspectroscopy (FS or F(z)):\n(i) For a BI monomer, FS were measured along a line intersecting two corner sites of the\nimidazole ring. Plotting the minimum value from each FS along this path generates an\nasymmetric F (d) curve with two peaks. These peaks correspond to nitrogen sites, with\nmin\nthe higher peak identifying the pyridine-like N site (Fig. 2d,e).\n(ii) Applying FS along the intermolecular bright line similarly differentiates NH and N sites\nwithin an H-bond (Fig. 2f-h). This method enables determination of H-bonding\ndirectionality in open chains (cf. Supplementary Fig. 6) and hexamer chirality (Fig. 3g).\n(4) The ability of FS to discriminate between N and NH sites originates primarily from atomic-\nscale electrostatic interactions between the CO tip and the two nitrogen species, which have\ndistinct electron density distributions.\n(i) To support this interpretation, DFT calculations of the Hartree (electrostatic) potential\nwere conducted for the -8.5°-L hexamer. The Hartree potential map reveals that the N site\nexhibits a more negative potential than the NH site (Fig. 3e). Since the terminal oxygen\natom of the CO tip carries a negative charge3,4, this results in a more repulsive interaction\nover the N sites.\n(ii) Using the PPAFM code, we simulated both AFM images and site-dependent ∆f(z) curves.\nThe characteristic ∆f (d) curves for six H-bonds and the ∆f(z) curves at the N and NH\nmin\nsites in the hexamer were successfully reproduced (Fig. 3h-j ), showing good agreement\nwith our experimental data (Fig. 3g).\nproduces an asymmetric, double-peaked F (d) curve, which can be used to determine the H-\nmin\nbonding directionality. We anticipate that this approach can be generalized to characterize H-bonds\nother than the N−H···N type.\nAction: Related discussion has been enhanced in the revised manuscript (page 6) and the simulated\nelectrostatic potential images as well as ∆f (d) curves for H-bonds in the hexamer have been\nmin\nincluded into the main text and SI (Fig. 3 and Supplementary Fig. 21).\nCommon question 2: The molecules might be slightly tilted on the surface, and the local\nadsorption environment could contribute to the difference in F . (Reviewers #1, #2 and #4)\nmin\nReply: We address this comment as follows:\n(1) DFT calculations based on experimentally determined adsorption registries reveal a flat-lying\ngeometry of the hexamers (Supplementary Figs. 11 and 20).\n(2) To exclude the possibility that variations in adsorption environments influence the measured\nF (d) curves for the BI hexamers, additional simulations were performed and compared with\nmin\nexperimental findings.\n(i) First, we examined the simulated and measured F (d) curves for the -8.5°-L and -8.5°-\nmin\nR hexamers (Supplementary Fig. 34). These configurations correspond respectively to\nstates 1 and 3 (Fig. 5e,f), and share identical BI backbone adsorption geometries but differ\nin H-bonding directionality. Remarkably, the two sets of F (d) curves for the cyclic H-\nmin\nbonds, simulated along clockwise paths, clearly reveal opposite H-bonding\ndirectionalities and closely match with experimental observations (Supplementary Fig.\n33), thus validating our line-FS approach.\n(ii) In addition, we modeled the -40°-L hexamer, whose experimentally determined\nadsorption geometry differs from that of the -8.5°-L hexamer (Supplementary Figs. 8).\nThe simulated F (d) curves again show a pronounced asymmetry associated with H-\nmin\nbonding directionality, consistent with experimental results (Supplementary Figs. 19 and\n21).\nIn summary, for all investigated BI hexamers, the observed asymmetry in the F (d) curves\nmin\nreflects the H-bonding directionality, whereas differences in adsorption geometry play only a\nminor role.\nAction: The SI now includes DFT-based PPAFM simulations of hexamers with varying chiralities\nand adsorption registries (Supplementary Figs. 21 and 34).\nCommon question 3: How does temperature affect proton ordering? Why chose simulation\ntemperature at 50 K? (Reviewers #2, #3 and #4)\n(1) The experimental evidence regarding the effect of temperature on proton delocalization\nremains unknown, as such experiments would require dedicated efforts beyond the scope of\nthis study.\n(2) Performing PIMD simulations at 50K represents a compromise between computational cost\nand the need to capture nuclear quantum effects (NQEs) at low temperature (LT). Conducting\nPIMD simulations at the DFT level at the experimental temperature (4.8 K) would demand\ntens of beads to achieve convergence. This would result in a computational cost beyond our\ncurrent capabilities.\nAlthough 50 K is higher than the experimental temperature, it still falls within a LT regime\nwhere thermal effects are negligible and proton delocalization remains dominated by quantum\nfluctuations. This approach has been validated in previous work, where PIMD simulations at\ntetramers on NaCl(001)5, consistent with experimental data obtained at 5 K6.\nTo further verify the reliability of our approach, we performed additional PIMD\nsimulations using LAMMPS7 and reaction force field8 at 50K and 5K, respectively. As shown\nin Supplementary Fig. 27, the resulting probability profiles at the two temperatures are\nqualitatively consistent, indicating that 50K lies within the deep tunneling regime. Therefore,\nusing 50K PIMD simulation at the DFT level to interpret the experimentally observed quantum\ndelocalization at 4.8K is reasonable.\nAction: PIMD simulations using LAMMPS package at 50K and 5K have been provided in the SI\n(Supplementary Fig. 27). The PIMD simulation of a hexamer on the Ag(111) surface, performed\nwith LAMMPS, is shown in Fig. 4d-f of the revised manuscript.\nCommon question 4: Whether an isotope effect would influence the findings? Performing\nexperiments with deuterated molecules could help discern conventional from delocalized\nhydrogen bonds and verify the role of quantum effects. (Reviewers #1 and #3)\nReply: To study the isotope effect, the deuterated BI (D-BI) molecules were synthesized (for\ndetails, see SI).\n(1) We first verify the deuteration by measuring the onset bias voltage for triggering the collective\nproton transfer in D-BI hexamers. The voltage threshold was found to be lower than that of H-\nBI hexamer (Supplementary Fig. 37), consistent with the reduced N-D stretching energy.\nLine-FS characterization of D-BI hexamer revealed both symmetric and asymmetric\nF (d) curves, which resemble those of H-BI hexamers (Supplementary Figs. 38 and 39). The\nmin\ncorresponding ∆F values for cyclic H-bonds show the same trend for D-BI and H-BI\nmin\nhexamers. Near F , AFM images for D-BI and H-BI hexamers display coexisting single- and\nmin\ndouble-lobe contrast.\n4\n(2) DFT-based PIMD simulations of a gas-phase D-BI yields a deuteron probability density\ndistribution and free-energy profile comparable to those of the H-BI dimer (Supplementary\nFig. 40), indicating deuteron delocalization in the H-bond.\n(3) Together, the experiments and simulations show that H/D substitution does not suppress\nquantum delocalization in BI hexamers. This is consistent with prior observations of deuteron\ndelocalization in Zundel cations (D O +) under both surface adsorbed9 and gas-phase10 UHV\n5 2\nconditions, as well as in low-barrier hydrogen bond between amino groups11.\nRuling out alternative origins\n(1) Tip effect: Repeating the measurements with a different qPlus sensor yields consistent results\n(data not shown), excluding tip-related artifacts.\n(2) Substrate effect: Comparable rotational-symmetry breaking in proton ordering was also\nobserved for BI hexamers on Au(111) (see Reply to Comment 1 from Reviewer #3),\nconfirming the generality of this phenomenon.\n(3) Geometric effect: We also checked the geometric effect by modeling a hexamer with\nintentionally varied H-bond lengths. After DFT relaxation, the H-bond lengths become\nequivalent (Supplementary Fig. 26).\nTo summarize, the rotational-symmetry breaking of proton ordering in cyclic H-bonds is intrinsic\nand reproducible across different qPlus sensors and substrates. The single- and double-lobed AFM\ncontrast can be altered by tip manipulation, yet this behavior is not reproduced by DFT simulations\nthat neglect NQEs. When NQEs are included, first-principles-based dynamic simulations can\nappropriately describe the proton delocalization observed in experiments, supporting NQEs as the\nmost natural and probable origin.\nAction:\n(1) The title has been revised to “Anomalous symmetry breaking in proton ordering of surface\nconfined cyclic hydrogen bonds” to emphasize the experimentally established findings, while\nreserving the NQEs attribution for the body of the paper.\n(2) Figure 1 has been updated: the scheme now shows that, for non-classical H-bonds, the proton\nis not equally shared between two nitrogen atoms, consistent with Fig. 4e,f.\n(3) Experimental line-FS data for D-BI hexamers have been added to the SI (Supplementary Figs.\n38 and 39), with related discussion in the “Discussion” section.\n(4) PIMD calculations for D-BI dimer have been added to the SI (Supplementary Fig. 40).\n“…The quality of the data is very good. However, the interpretation of the results and the claimes\nmade are not very convincing. I do not think that there is enough experimental evidence for the\ndistinction between classical and non-classical hydrogen bonds. It would be great, but is not\nconvincingly shown. Several questions remain open and need to be considered. I recommend that\nthe paper is reconsidered after mayor revision.”\nComment 0: Overall, the data is of high quality. However, I find the assignment of classical\nhydrogen bonds (N−H···N) and non-classical hydrogen bonds with quantum proton delocalization\n(N···H···N) not convincing. I also see a potential problem with hexamers, that show both type of\nhydrogen bonds. If molecules have a classical hydrogen bond on one side, and thus either have a\nhydrogen or no hydrogen bonded to that N atom, it requires the other N atom of that molecule to\nhave no or a hydrogen, respectively, and thus also a classical hydrogen bond, or they would\nbecome radicals, right? I can imagine a structure with all non-classical hydrogen bonds. However,\nin the mixed hexamer, with both types of bonds, I see a conceptual problem. This needs to be\ndiscussed in more detail, or maybe reconsidered.\nReply: We appreciate the reviewer’s positive assessment of our data quality. As this comment is\nsimilar to Comment 4, we provide a detailed response in the Reply to Comment 4 from Reviewer\n#1.\nComment 1. I think it is also very instructive to consider the closed shell Kekule structure of\nbenzimidazole. In it, the two bonds that connect the carbon atom at the 5 membered ring to N and\nNH site are different. A formal double bond to the N and a formal single bond to the NH group.\nThe resulting bond-order related contrast with AFM is in line with the AFM image of the chain\nshown in Fig. 2f. The formal double bonds appear shorter compared to the single bonds. Moreover,\nthe related difference in electron density above the two N atoms might contribute to the contrast\nmechanism on N and NH sites, as observed in the monomer. The bond order of N-C bonds in the\nhexamers could be analysed.\nReply: Following the reviewer’s suggestion, we analyzed the apparent bond lengths of the single\n(C−N) and double (C=N) bonds in the imidazole rings of BI molecules in both open chains and\nhexamers. The results are summarized below.\n(1) Open chains: Since the H-bond directionality in open chains can be determined using FS\nmeasurements, the formal single and double bond positions in the imidazole ring can be\nidentified. To quantitatively assess the bond-order contrast, we defined a ratio γ as (𝐿 −\n𝐶−𝑁\n𝐿 )/ 𝐿 , where 𝐿 and 𝐿 represent the apparent lengths of the C−N and C=N\n𝐶=𝑁 𝐶−𝑁 𝐶−𝑁 𝐶=𝑁\nbonds, respectively. Supplementary Fig. 7 summarizes the γ values from two 9-membered\nchains, where the predominance of positive γ values indicates that 𝐿 is consistently larger\n𝐶−𝑁\nthan 𝐿 .\n𝐶=𝑁\n6\n(2) Hexamers: Similarly, FS measurements enabled determination of H-bond directionality in\nhexamers and thus the identification of C−N and C=N bond positions within imidazole rings.\n(i) In contrast to the linear chains, negative γ values were found for the H-BI (or D-BI)\nhexamers (Supplementary Figs. 12 and 13), implying that the apparent length 𝐿\n𝐶=𝑁\noccasionally exceeds 𝐿 .\n𝐶−𝑁\n(ii) However, no systematic trend in γ values for hexamers was identified. First, within a\nsingle hexamer, the γ values display both positive and negative signs. This makes it\ndifficult to determine the chirality (R or L) of a hexamer solely from bond lengths.\nMoreover, for the same type of hexamers (not necessarily the identical one), repeated\nmeasurements often yield variations in bond lengths or in the sign of γ (Supplementary\nFig. 12i,k or j,l), likely due to slight differences in tip termination. In contrast, force\nspectroscopy measurements appear highly robust against such tip-dependent variations.\nIn the -8.5°-L-1 hexamer (Supplementary Fig. 12a,b), FS measurements reveal more\nsymmetric F (d) curves for H-bonds involving molecule L3, L4 and L5 (highlighted in\nmin\ngreen in Supplementary Fig. 12i), suggesting proton delocalization. One would therefore\nexpect smaller γ values for these molecules compared to L1, L2 and L6. However, such\na trend is absent in both -8.5°-L-1 and in -8.5°-L-2 (Supplementary Fig. 12i,k). Similarly,\nfor the +8.5°-R-1,2 hexamers, despite FS confirming delocalized protons between\nmolecules R2 to R5, neither the sign nor magnitude of γ exhibited consistent behavior\n(highlighted in green in Supplementary Fig. 12j,l).\n(iii) Applying the same analysis to D-BI hexamers reveal no systematic trend in γ. However,\nsome molecules in the hexamer show larger |γ| (Supplementary Fig. 13).\nTherefore, we conclude that compared to linear chains, the bond-order contrast in hexamers is\nmore irregular and thus does not serve as a good indicator of proton delocalization in this system.\nAction: The related analysis bond-order associated apparent bond length contrast for both linear\nchains and hexamers have been provided in the SI (Supplementary Figs. 7,12 and 13). Related\ndiscussion has been included into the main text (page 5,6).\nComment 2: In the hexamer the differences in bond-order contrast are less compared to the open\nchain. This should be discussed. Might it imply that possibly in all hexamer structures all the H\nbonds are delocalized?\nReply: For the bond-order contrast analysis of hexamers, please see our Reply to Comment 1\nfrom Reviewer #1. Line-FS characterization on hexamers resolve both classical and non-classical\nH-bonds (details in our Reply to Comment 4 from Reviewer #1).\nComment 3: In a non-classical hydrogen bond with quantum proton delocalization (N···H···N),\nwould it need to be symmetric? Could not the proton be delocalized, but still be (in average) closer\nto one of the two N, allowing some asymmetry of the N···H···N bond?\nnamely a hexamer adsorbed on Ag(111). The results reveal that the proton delocalization in\nshortened H-bonds is strongly asymmetric (Fig. 4e,f). Details are provided in Reply to Comment\n4 from Reviewer #1.\nComment 4. A structure that has in part classical hydrogen bonds (N−H···N) and in part non-\nclassical hydrogen bonds with quantum proton delocalization (N···H···N) as claimed in states 1\nand 4 in Fig. 1 seems problematic to me. The molecules connecting non-classical and classical H\nbonded segments have a “wrong” count in hydrogen atoms. That is having one classical hydrogen\nbond in the ring should determine all other classical hydrogen bonds in the ring. It seems to me that non-classical\nhydrogen bonds with quantum proton delocalization (N···H···N) should only be possible for either all bonds in the hexamer,\nor none. Please discuss in detail how a molecule that on one side has a classical H bond, can have a delocalize one on its\nother side.\nReply: We rationalize the coexistence of classical and non-classical H-bonds in the hexamer as\nfollows:\n(1) We performed LAMMPS-based PIMD (see methods) for the full system: a hexamer adsorbed\non the substrate. To reproduce the experimentally observed symmetry breaking in proton\nordering, we introduced a geometrical distortion resulting in varying H-bond lengths. During\nthe PIMD simulations, only the hydrogen atoms were permitted to move.\n(i) Beginning with the DFT-optimized structure (Fig. 3b), which features threefold rotational\nsymmetry and nearly equal H-bond lengths, we induced asymmetry by rotating molecules\nL2, L3, L5, and L6 toward L4 around the hexamer center at a fixed radius, sweeping an\nangle of |4.8°| (Fig. 4d; Supplementary Fig. 27a). This adjustment reduced N···N\nseparations for L3-L4 and L4-L5 to ≈ 2.4 Å, maintained L2-L3 and L5-L6 near 2.8 Å,\nand elongated L1-L2 and L6-L1 to ≈ 3.18 Å.\n(ii) Proton delocalization in H-bonds strongly depends on bond length12, and occurs only in\nthe shortened H-bonds L3-L4 and L4-L5 (Fig. 4d). In contrast, for N···N ≳ 2.8 Å, the\nproton remains localized at the donor atom (Supplementary Fig. 27b,c).\n(iii) A key feature of the shortened H-bonds is the pronounced asymmetry in both the proton-\ndensity distributions and the associated free-energy profiles (Fig. 4e,f). This indicates that\nthe proton primarily localizes near the donor site, with substantially lower probability at\nthe acceptor site of the neighboring molecule. This reaction is dominated by quantum\ntunneling, thereby reducing the effective free-energy barrier to a few meV.\n(iv) We further compare the proton delocalization in the gas-phase BI dimer with that in the\nadsorbed hexamer. In the gas-phase dimer, removing one proton results in two mirror-\nsymmetric configurations that interconvert via proton tunneling, giving rise to a\nsymmetric probability distribution and a double-well potential (Fig. 4a,b).\nIn comparison, in the hexamer the molecule L5 retains a fully localized proton on its\nleft side, whereas L3 is proton-deficient on its right side. As a result, proton tunneling\n8\nacross H-bonds L4-L5 and L3-L4 transiently generates a benzimidazolium cation (L5)\nand a radical benzimidazole species (L3). This asymmetry leads to unequal proton density\ndistributions and asymmetric free energy profiles, thus affecting the degree of proton\ndelocalization (Fig. 4e,f).\n(2) Imidazole is amphoteric and undergoes autoprotolysis to a much higher degree than water\nleading to high proton conductivity13. Notably, imidazolium cation (i.e., the conjugated acid\nof imidazole) has a pK of approximately 7.0 in aqueous solution14. This means that in neutral\na\nwater, roughly 50% of the imidazole molecules are protonated. Moreover, gas-phase\ncomputational studies reveal a proton affinity of ~941 kJ/mol at N3, significantly higher than\nat other positions (N1 or the ring carbons)15. Protonation at N3 is therefore both kinetically\naccessible and thermodynamically favored due to resonance stabilization and aromaticity\npreservation.\nIn summary, although surface-stabilized protonated or deprotonated BI species were not\ndirectly observed in the intact hexamer, our LAMMPS-based PIMD simulations indicate that\ntransiently formed protonated (e.g., L4, L5) and deprotonated (L3) BI species coexist due to\nquantum proton delocalization in shortened H-bonds, albeit with lower delocalization probabilities.\nAccording to our DFT-based PIMD simulations, the quantum proton tunneling occurs on a\ntimescale below one picosecond (cf. Fig. 4c). Consequently, time-averaged FS measurements\ncannot directly capture these transient benzimidazolium cation and deprotonated BI species in the\nhexamers.\nAction: PIMD simulations of the hexamer on the Ag(111) substrate using the LAMMPS code\nhave been added to Fig. 4d-f, along with corresponding text in the revised manuscript (page 8).\nThe Fig. 1 has also been revised to place the proton off-center of in the non-classical H-bond,\nreflecting its asymmetric probability distribution.\nComment 5. Related to the previous point: In theory the hydrogen bonds with quantum proton\ndelocalization (N···H···N) is only obtained if one hydrogen is removed from one molecule (Fig.\n4b). It seems this “trick” is needed because of the problem described in the previous point. What\nis the charge state of the structure in Fig. 4b? Is it a neutral doublet, or charged?\nReply: In the gas-phase dimer model we removed one hydrogen atom and added an extra electron,\ncreating a negatively charged system in which two proton-transfer states are mirror-symmetric,\nyielding a symmetric double-well free-energy profile. By contrast, in the complete system of\nadsorbed hexamer, proton delocalization arises naturally once the H-bond length is sufficiently\nshort (Supplementary Fig. 27), without a charge adjustment, even though those calculations were\nperformed using a classical force field.\nAction: Details for the dimer setup for DFT-based PIMD calculations are included in the SI (page\n6).\npossible explanation for the observed switching and asymmetric contrast in the hexamers (and\nrelate to the different 8.5 and 40 deg tilted structures).\nReply: Removing a single hydrogen atom from the hexamer significantly disrupts its balanced and\nsymmetric ring structure.\nComment 7. The authors make not clear how they distinguish between classical hydrogen bonds\n(N−H···N) and non-classical hydrogen bonds with quantum proton delocalization (N···H···N).\nMore specific: How do the authors come to the specific value of 1 pN, that they use to discern\nclassical and non-classical N bonds? Why this specific value? Possibly, (N···H···N) bonds\nmight be asymmetric. Moreover, the molecules might be tilted slightly on the surface or their local\nadsorption environment might be different possibly contributing to the Fmin difference, that is not\nrelated to the nature of the hydrogen bond. The value of 1 pN to discern seems chosen arbitrary.\nReply: We address the raised questions as follows:\n(1) For better distinguishing between classical and non-classical H-bonds in BI hexamers, we\nsimulated AFM images and ∆f (d) curves using DFT-optimized structures and DFT-\nmin\ncalculated Hartree potential, combined with PPAFM package (see Reply to Common\nquestion 1).\n(i) AFM images: In simulation, all six BI molecules exhibit a bright single-lobe feature on\nthe right side (Fig. 3f). In contrast, the experimental image displays both single- and\ndouble-lobe features within the hexamer (Fig. 3c).\n(ii) ∆f (d) curves: For all six H-bonds, the simulated curves consistently display two peaks\nmin\nwith unequal heights (Fig. 3h). Experimentally, specific bonds (L3-L4 and L4-L5) show\nnearly symmetric M-shaped curves.\nThese discrepancies provide a criterion: H-bonds with nearly symmetric M-shaped F (d)\nmin\ncurve and double-lobe AFM contrast near F are assigned as non-classical.\nmin\n(2) We define |∆F | = |F (NH) − F (N)|. Reference values are as follows:\nmin min min\nIsolated BI monomer: |∆F | ≈ 7 pN (Fig. 2e);\nmin\nOpen chain: |∆F | ≈ 4 pN (Fig. 2h and Supplementary Fig. 6);\nmin\nSome H-bonds in hexamers: |∆F | ≳ 4 pN (Supplementary Fig. 15i,j).\nmin\nWe adopt |∆F | ≲ 1 pN is a conservative upper bound for assigning a non-classical H-bond:\nmin\nit lies well below all classical references and within our FS resolution. A looser cutoff of |∆F |\nmin\n≲ 2 pN would leave our assignments unchanged except that L2-L3 in −8.5˚-L would also be\nclassified as non-classical (Supplementary Fig. 15i), consistent with its faint double-lobe AFM\ncontrast (Supplementary Fig. 15e).\n(3) For possible molecular tilt and adsorption contribution to asymmetric F (d), please refer to\nmin\nReply to Comment question 2.\nAction: We incorporate the text describing our criterion for distinguishing classical and non-\nclassical H-bonds in the SI (page 4).\nwith a hcp site (Fig. 3d and Supplementary Fig. 7). This indicates that the mirror symmetry for\nfreestanding 6-membered enantiomorphs does not hold when considering their adsorption\nregistries (Fig. 3c,d).” I do not understand the argument. If the center would not be the same site,\nI would agree that the structures cannot be mirror symmetric. I do not see why both structures\nbeing centred around hcp sites makes the authors conclude they are not mirror symmetric. In fact\nthey do look mirror symmetric to me. (Flip at horizontal line in Fig. S7).\nReply: The reviewer is correct. Flipping the −8.5°-L-hexamer horizontally together with its\nsubstrate, generates the +8.5°-R-hexamer without changing the adsorption registry.\nAction: The sentence has been revised to “This indicates that the −8.5°-L-hexamer and the +8.5°-\nR-hexamer are mirror symmetric with respect to the 〈1̅1̅2〉 directions”.\nComment 9. Related to the previous point: Is the energy difference calculated for the enantiomers\nin Fig. S8 significant? Or is it within the noise (energy cutoff criterium) of the calculation?\nReply: The total absorption energy difference of the two enantiomers is 8 meV, which is an\nextremely small value difference that is beyond the accuracy of DFT calculations and falls within\nthe typical numerical noise of DFT calculations.\nAction: We have added a sentence in the SI (Supplementary Fig. 11) to clarify this point.\nComment 10. The energies for adsorption under +/-8.5 deg (Fig. S8) could be compared with\nsimilar calculations for the local adsorption minima expected at +/-40 deg.\nReply: We respond to this suggestion with the following clarifications:\n(1) We calculated the adsorption energy (E ) for ±40° hexamers centered at an hcp hollow site\nad\nand obtained E ≈ -5.72 eV. This value is ~0.15 eV lower than that of the ±8.5° hexamers (E\nad ad\n≈ -5.57 eV; Supplementary Fig. 20).\n(2) We also computed the cyclic H-bonds binding energy (E ) of the ±8.5° and ±40° hexamers and\nb\nobtained similar values, E ≈ -3.1 eV.\nb\n(3) Thus, DFT calculations indicate that ±40° hexamers are energetically more favorable on\nAg(111) compared ±8.5° hexamers. However, in experiments, ±8.5° hexamers are more\n11\nprevalent, implying a potential contribution from NQEs that are absent in DFT energy\ncalculations (Supplementary Fig. 20). Supporting this interpretation, DFT relaxations based on\nthe experimentally determined adsorption registry show that in ±8.5° hexamers, pyridine-like\nN atoms sit close to hollow site, a geometry that facilitates proton tunneling (Supplementary\nFig. 35).\nAction: The corresponding DFT calculation and related discussion have been included into the SI\n(Supplementary Fig. 20).\nComment 11. The authors write: “Thus, our methodology demonstrates that H-bond directionality\ncan be unambiguously determined via line-FS measurements, offering a clear resolution to the\nlongstanding debate over whether the bright intermolecular contrast represents an H-bond or an\nartifact due to CO-tip flexibility [31,36].” First of all: what is the resolution to the debate?\n(I fear the authors might mean that the contrast is due to the hydrogen bond). What is the reasoning? In\nthis context, please consider that you also state: “This feature was observed in the Fmin(d) curves\nmeasured for a single BI molecule and linear chain (Fig. 2e,h), where Fmin(NH) is lower than\nFmin(N).” Note, that you observed the feature also in a BI monomer, that is in the absence of any\nhydrogen bonds.\nReply: We provide our response to this key concern in the Reply to Common question 1.\nComment 12. Is the adsorption site of the structures shown and switched in Fig. 5 identical? That\nis, are the transition states 2 and 3 also centered around hcp sites?\nReply and Action: Transition states 2 and 3 also center around the hcp sites. Corresponding\nexperimental data have been added into the SI (see Supplementary Fig. 29).\nComment 13. Related to previous point: Can the 40 deg oriented structures also be switched or\nmanipulated to become other structures, e.g. the +/-8.5 deg structures?\nReply and Action: Yes, the +40° hexamer can be obtained via tip manipulating of the −8.5°\nhexamer. Additional experimental data has been included into the SI (Supplementary Fig. 28d,e).\nComment 14. If orbital resonances could be accessed by STM, might that provide additional\nexperimental data on the type of hydrogen bond?\nReply: To address this point, we performed additional scanning tunneling spectroscopy (STS or\ndI/dV) measurements on a single hexamer (Supplementary Fig. 23). (i) Overview STS spectra,\nrecorded over a ±2.2 V bias range, were taken on individual BI molecules within the hexamer\n(Supplementary Fig. 23a). Near the Fermi level, the spectra are featureless, primarily reflecting\nthe Ag(111) surface state. A pronounced increase in the dI/dV signal is observed above +2V, which\ncan be ascribed to the detection of molecular orbital resonances located far away from the Fermi\n12\nenergy. (ii) For further inspection, we performed dI/dV mapping at three selected bias voltages.\nAs shown in Supplementary Fig. 23b-d, the resulting LDOS maps are spatially homogenous and\ndo not reveal localized features associated with H-bonds. Therefore, we conclude that nc-AFM is\na more suitable method, which is sensitive to the localized electron density.\nAction: dI/dV spectra on BI molecules in hexamer have been included in the SI (Supplementary\nFig. 23).\nComment 15. Optional: The paper could very likely be greatly improved by performing\nexperiments with deuterated molecules. Effects of proton tunneling could be revealed very\nconvincingly by the comparison of proton vs. deuteron tunneling. It could help discerning\nconventional from delocalized hydrogen bonds. In the present manuscript the assignment and\ndistinction of these two bonds is not clear.\nReply: We provide a detailed response to this key concern in the Reply to Common question 4.\nPoint-by-point response to the Reviewer #2’s comments\n“…This work is highly relevant to the study of quantum effects in hydrogen bonds and\nmolecular self-assembly on surfaces. The manuscript provides an important insight into\nquantum proton ordering and hydrogen bonding at the atomic scale. However, to ensure the\nvalidity of the conclusions, the authors should address major concerns. Additionally, minor\ncorrections and refinements in figures and terminology are necessary.”\nComment 1-i. The primary issue with the manuscript is the interpretation of hydrogen bond\ndirectionality based on force spectroscopy. As the authors are likely aware, there is ongoing debate\nin the community regarding whether the contrast observed above hydrogen bonds originates from\nthe bonds themselves or other effects. Specifically, it remains unclear whether CO-tip AFM\ntruly images hydrogen bonds or if the observed contrast stems from other intermolecular interactions\n(see, e.g., S. K. Hämäläinen et al., Phys. Rev. Lett. 113, 186102 (2014)). Since this assumption is\nfundamental to the manuscript's conclusions, the lack of supporting data beyond point force\nspectroscopy is surprising. A straightforward way to validate these findings would be to provide\na series of BR-AFM images at different tip-sample distances to corroborate the spectroscopy results.\nReply: We address the questions as follows:\n(1) For a detailed explanation of how force spectroscopy reveals H-bond directionality, please\nrefer to Reply to Comment questions 1.\n(2) As requested, we carried out height-dependent AFM mapping in greater detail (Supplementary\nFig. 18). By analyzing AFM slices at different heights, we extracted the force spectrum, Δf(z),\nat selected sites (x,y). Along the apparent H-bond, we traced a line to obtain the Δf (d) profile\nmin\n13\n(Supplementary Fig. 18b). The curve exhibits an M-shaped feature, with the values at the N\nsite consistently higher than those at the NH site.\nFor several representative heights corresponding to the force spectra, the associated AFM\nimages were obtained (Supplementary Fig. 18f). At the height of the force minimum (most\npronounced at z₂), the N site appears brighter, corresponding to the single lobe in the AFM\nimage. This analysis thus links the experimental force spectra with the AFM images,\ndemonstrating that the lobe observed in the image originates from the force contrast between\nthe N and NH sites. Furthermore, the experimental and simulated force spectra and AFM\nimages show excellent agreement (Supplementary Fig. 18b-g), thus confirming the validity of\nour approach.\nAction: Height-dependent AFM images and related analysis have been added into the SI\n(Supplementary Fig. 18). The paper by Hämäläinen et al. has been cited in the revised manuscript.\nComment 1-ii. Additionally, the Probe Particle Model is a well-established method for simulating\nBR-AFM images with high accuracy. Given that this code is open-source, the authors could\nperform their own simulations to support their interpretations. Interestingly, my own simulations\nof BI in the gas phase did not reveal significant differences between the N and N-H sites,\nsuggesting that the observed contrast might be a topographical effect rather than a direct measure of\nhydrogen bonding. Further scrutiny of the provided images suggests that the molecule may not\nbe lying flat on the surface, which could influence the total force measurements.\nReply: We thank the reviewer for this suggestion. For comprehensive discussion on PPAFM\nsimulation and the effect of molecular tilt on the force measurements, please refer to Reply to\nCommon questions 1 and 2.\nComment 1-iii. Another noteworthy observation is the discrepancy in the minimum force values\nextracted from Figures 2e and 2h, which range from ~10 pN to ~40 pN. The authors should clarify\nthe physical origin of this variation—could it be due to incomplete background signal removal?\nReply: This discrepancy is not due to incomplete background subtraction. The detailed procedure\nof deriving 𝐹(𝑧) is described in the Supplementary Fig. 16: we first measured the force spectra on\nthe bare Ag(111) surface up to a sufficiently large tip–sample distance where Δ𝑓 becomes nearly\nconstant. Using the same parameters, we then measured the spectra on the molecule. At large\ndistances, the spectra acquired on the surface and on the molecule overlap, allowing background\nsubtraction. The residual frequency shift was subsequently converted to force using the standard\nSader–Jarvis method16.\nTo rule out artifacts, we repeated FS measurements on both monomers and hexamers.\nConsistently, the minimum force at the N/NH sites of a single molecule is ~ −20 pN, whereas\nwithin hexamers it reaches ~ −40 pN. This systematic difference is intrinsic and stems from distinct\nchemical environments: an isolated molecule lacks intermolecular interactions, while in linear\nchains or hexamers each molecule engages in multiple H-bonds. This alters the local charge\ndistribution and manifests directly in the measured forces.\nto-noise ratio, and these new data have replaced the previous measurements in both the main text\n(Fig. 2e) and the SI (Supplementary Fig. 4).\nComment 1-iv. Additionally, why should the N-H site be more attractive to the CO tip than the N\nsite? This could be better understood by analyzing the electrostatic potential (ESP) of the molecule\nor performing density functional theory (DFT) calculations of the force spectroscopy. The authors\ndiscuss the role of adsorption registry in proton ordering, yet the specific electronic interactions\nbetween the substrate and the molecules remain unclear. Does the Ag(111) surface actively\nstabilize particular proton configurations? Are charge transfer effects involved in influencing\nhydrogen bonding? Additional discussion or computational analysis would strengthen these claims.\nReply: We respond to this comment as follows:\n(1) We have computed ESP of the full system at the DFT level. A detailed discussion can be found\nin Reply (point 4) to Common question 1.\n(2) To clarify how adsorption geometry governs proton ordering in the cyclic H-bonds, we\ncalculated electron density difference (Δρ) for the -8.5°-L (state 1) and +8.5°-L (state 2)\nhexamers on Ag(111) (Supplementary Fig. 35).\nIn the -8.5°-L (state 1) hexamer, which favors the proton delocalization, all pyridine-like\nN atoms occupy hollow sites, whereas in the +8.5°-L (state 2) hexamer, three pyridine-like N\natoms occupy top sites (Supplementary Fig. 32c). Cross-sectional Δρ maps taken along the H-\nbonds show a much larger charge redistribution beneath proton-acceptor nitrogen atom (N’) at\ntop sites (Fig. S35d) than beneath its hollow site counterpart (Fig. S35c). This enhanced\nelectronic hybridization at top sites markedly tilts the nitrogen σ-orbital (lone pair) out of the\nmolecular plane (Fig. S35d), thereby suppressing proton tunneling in state 2 (Supplementary\nFig. 32g).\nAction: Related discussion can be found in the main text (page 10, the paragraph before\nDiscussion), with corresponding figure in Supplementary Fig. 35.\nComment 2. The manuscript suggests that proton transfer in BI hexamers is primarily induced by\ninelastic electron excitation from the STM tip. However, alternative mechanisms, such as\nthermally activated or substrate-mediated processes, cannot be ruled out. Have the authors\nconducted control experiments at varying temperatures or tip biases to exclude these possibilities?\nCan the authors also provide a mechanism explaining how inelastic electron tunneling induces\nproton transfer?\nReply: We address the raised questions as follows:\n(1) Additional mechanisms for proton transfer within the hexamer cannot be ruled out. Thermal\nactivation at 77 K may contribute. As demonstrated in Supplementary Movie1, tip-induced\ndisturbances during scanning readily trigger hexamer rotation between two mirror-symmetric\n15\nadsorption configurations, indicating interconversion between enantiomorphic structures via\ncollective proton transfer.\n(2) To explore the mechanism of inelastic electron tunneling-induced proton transfer, we\ncalculated vibrational modes for a BI monomer absorbed on Ag(111). The NH stretching mode\nwas found at ~389 meV (Supplementary Fig. 31), closely matching the experimentally\nobserved threshold voltage. This indicates that collective proton transfer in the hexamer is\ntriggered when the tunneling electron energy coincides with the NH stretching energy. A\nsimilar mechanism has been reported in intramolecular hydrogen atom transfer reactions of\nsurface-adsorbed single porphycene molecules17.\nAction: The mechanism is provided in the main text (page 9), with additional calculation\n(Supplementary Fig. 31) and a reference (Ref. 47) to support inelastic-electron-tunneling triggered\nproton transfer.\nComment 3. Since all experiments are conducted at 5 K, it is unclear how temperature affects\nproton delocalization and tunneling. Would quantum proton ordering persist at higher\ntemperatures, or would classical hydrogen bonding become dominant? A discussion on this aspect\nis necessary.\nReply: We computationally investigated the effect of temperature on proton delocalization by\ncomparing the proton motions in hexamer on the surface from PIMD simulations with classical\nforce field at 5K and 50K. A detailed discussion of these results is provided in Reply to Common\nquestion 3. Conducting temperature-dependent experiments to probe proton ordering requires\nsubstantial additional effort and it lies beyond the scope of this study.\nAction: Comparative simulations at 5 and 50 K have been included into the SI (Supplementary\nFig. 27).\nComment 4. Figures 3 and 4 require improved annotations to clearly indicate key measurements,\nparticularly the extracted force-distance parameters.\nReply: We have updated Figs. 3 and 4. In the revision we define two tip-sample distance variables:\nz’, reference to the STM setpoint (used for imaging heights), and z, the displacement used in F(z)\nwith z = 0 referenced to the BR-AFM imaging height.\nAction: Figure captions in the main text have been revised accordingly, and the SI provides a more\ndetailed definition of z’ and z (page 5).\nComment 5. The phrase “broken rotational symmetry” appears frequently but is not well defined.\nThe authors should clarify its precise meaning and implications.\nUsing the experimentally determined adsorption registries (Supplementary Fig. 10),\nDFT\ncalculations for the ∓8.5°-L/R and ∓40°-L/R hexamers yield AFM contrast at F with nearly\nmin\nsixfold rotational symmetry, a flower-like pattern with six single lobes (Fig. 3f and Supplementary\nFig. 21c). By contrast, experimental AFM images acquired near F do not exhibit this symmetry\nmin\nbecause single- and double-lobe contrasts coexist (Fig. 3c). This broken symmetry in proton order\nreveals two types of H-bonds in the ground-state hexamers: classical and non-classical.\nAction: We have added a clearer definition of rotational-symmetry breaking to the main text (page\n7).\nComment 6. The introduction should provide a more detailed comparison with previous studies\non quantum proton delocalization in similar molecular systems.\nReply and Action: Previous studies on quantum proton delocalization in interfacial molecular\nsystems has been provided in the introduction (page 3).\nComment 7. Typographical and Formatting Issues:\no Page 3, line 7: “proton delocalisation” → “proton delocalization” (for consistency with\nAmerican English).\no Page 7, line 23: “image” → “images”; “Large” should be lowercase.\no Page 7, lines 24-25: The word “products” appears redundantly.\no Figure 3F appears misaligned compared to Figure 3B—this should be corrected.\no Figure 5: The red arrow is not explained in the caption.\no Supporting Information: Some materials (e.g., Scheme 1, Figures S1 and S2) are not referenced\nin the main text. Their relevance should be clarified.\nReply: We thank the reviewer for the careful reading.\nAction:\no We adopted American English consistently throughout the manuscript.\no Page 7, line23/lines 24-25: We have changed the “image” to “images” in the caption of Fig. 2.\nFor the latter case, we carefully checked the manuscript and could not locate the places the\nreviewer referred to. Accordingly, no changes were made.\no Figure 3: we verified the layout; Supplementary Fig. 13b and 13f are properly aligned.\no Figure 5: The arrows are now explained in the caption.\no Figures 1, S2 and S3 are cited in the main text.\nPoint-by-point response to the Reviewer #3’s comments\n“…A major strength lies in the ability to manipulate and switch between quantum and classical\nproton-ordering states using an AFM tip, offering insights into proton transport and potential\napplications in quantum materials. In my view, this research makes a significant contribution\nto the understanding of molecular-scale quantum effects and merits publication, as it introduces\nboth novel insights and technical advancements in the field. However, I have several questions\nand comments that should be addressed before acceptance, which I list below.”\nComment 1. The study is restricted to benzimidazole hexamers on Ag(111), which may not fully\nrepresent biological or bulk hydrogen-bonded systems. Surface interactions with Ag(111) may\ninfluence proton behavior in ways that do not generalize to solution-based or gas-phase\nenvironments. I would appreciate a comment from authors whether the observed quantum proton\nordering extends beyond the Ag(111) surface? Would the same proton ordering be observed on\nweakly vs. strongly interacting surfaces? Also, would this behavior persist in a freestanding\nbenzimidazole hexamer or in bulk crystalline environments? In any case, I believe a short\njustification on the choice of Ag(111) surface would be appreciated by the reader.\nReply: To explore whether proton delocalization also occurs on other surfaces, we carried out\nadditional experiments on BI molecules adsorbed on Au(111), Cu(111) and bilayer\nNaCl(001)/Ag(111).\n(1) On Cu(111), BI molecules bind strongly to the surface and no regular assemblies were\nobserved (Supplementary Fig. 36a,b). We anticipated weaker adsorption on NaCl(001).\nHowever, BI molecules show strong molecule-substrate interaction, evidenced by tilted\nadsorption in molecular chains (data not shown). Consequently, no hexamers form on\nNaCl(001)/Ag(111) (Supplementary Fig. 36c,d).\n(2) On Au(111), a small number of BI hexamers were found coexisting with linear chains\n(Supplementary Fig. 36e,f). We measured line-FS across six apparent H-bonds, and the\nresulting F (d) curves are shown in Fig. S36i. A comparable proton ordering, featuring the\nmin\ncoexistence of classical and non-classical H-bonds, was also observed on this more inert\ncoinage metal surface (Supplementary Fig. 36h,i), supporting the generality of the proton\ndelocalization in surface-confined cyclic H-bonds.\n(3) (i) In the revision, we carried out PIMD calculations with classical force field on a constrained\nhexamer/Ag(111) system (Supplementary Fig. 27). These results suggest that for a\nfreestanding hexamer, proton delocalization is still expected if the in-plane constraint and the\nshortened H-bonds are maintained. (ii) Previous studies of imidazole crystals have shown that\nH-bonded chains can facilitate collective proton transfer, enabling room-temperature\nferroelectricity18. However, whether NQEs contribute to such collective proton dynamics in\nbulk crystal remain elusive. (iii) Finally, imidazole rings an important structural unit in\nbiological systems, notably in enzymatic catalysis, where it stabilizes enzyme-intermediate\ncomplexes via the formation of low-barrier hydrogen bond (LBHB)19,20, in which NQEs play\na predominant role.\nAction: We added proton ordering characterization of BI hexamer on Au(111) to the SI\n(Supplementary Fig. 36), with related discussion implemented in Discussion (page10).\nbehavior, raising questions about the true ground-state configuration in the absence of an external\nprobe. Could authors comment on that?\nReply: We address the question through the following points:\n(1) The rotation of the hexamer requires a metal tip, whereas CO-terminated tips are utilized for\nAFM characterization. It is known that nc-AFM using a qPlus sensor functionalized with a\nCO molecule is suitable for imaging weakly bonded systems21,22. For example, a combined\nexperimental and theoretical study of water clusters on the NaCl(001) surface demonstrated\nthat the very weak tip-sample interaction allows precise structural determination with minimal\ndisturbance4.\n(2) The force exerted on the hexamers are on the order of tens of pN, which can be considered as\nminute. Moreover, our criteria for distinguishing between classical and non-classical bonds\nrely on the minimum value of the force curve (F ), which corresponds to the attractive tip-\nmin\nsample interaction region and thus implies a large-tip sample distance.\nTherefore, we conclude that the probed the proton ordering in the hexamer corresponds to a\nground state configuration.\nAction: A clarifying sentence has been added to the main text (page 7) with new reference (Ref.\n41), stating that the proton order is unperturbed by a CO-tip.\nComment 3. How does the temperature as well as external electric field affect proton ordering?\nCould an external field stabilize or disrupt the quantum proton sites? In case of temperature, I\nwould appreciate a comment why temperature of 50K was chosen for simulations.\nReply: We respond to raised questions as follows:\n(1) For temperature effects and our choice of 50K for DFT-based PIMD simulation, please refer\nto our detailed response in Reply to Common question 3.\n(2) To test whether an external electric field perturbs proton ordering, we acquired proton-order\nsensitive imaging at low bias. No discernible change was observed (Supplementary Fig. 22),\nthus indicating that weak electric fields have no impact on the proton ordering.\nAction: Additional experimental data have been included into the SI (Supplementary Fig. 22). We\nalso added related discussion into the main text (page 7).\nComment 4. Have authors consider using other techniques to further confirm that the observation\nis truly quantum mechanical behavior? For example, inelastic neutron scattering or high-resolution\nvibrational spectroscopy could confirm the quantum vs. classical nature of the discussed proton\nbehavior.\nspectroscopy could provide complementary evidence for quantum proton behavior. However,\nthese measurements require dedicated instrumentation and expertise that lie beyond the scope of\npresent study. Although imidazole-related compounds had been investigated by the\naforementioned techniques23-26, NQEs in hydrogen-bonded imidazole systems have not yet been\nreported. We therefore hope that our results will simulate further work in this direction using\ndifferent experimental approaches.\nComment 5. I also wonder whether an isotope effect would influence the findings? Would for\nexample deuteration lead to observable shifts in proton ordering, verifying the role of quantum\neffects? Could authors comment on that?\nReply: To address this key concern, we have performed additional experiments on deuterated BI\nmolecules. Please refer to Reply to Common question 4 for our detailed response.\nAction: The relevant discussion is now included in the Discussion section, supported by new data\nand PIMD calculations (Supplementary Figs. 38-40).\nComment 6. On a smaller note. Perhaps I am exaggerating, but in Fig. 4b and 4c, I notice that the\nPIMD simulation results appear to show two maxima/minima on the acceptor side. I am curious\nwhether this is simply an artifact of the simulation or if it could indicate additional effects that\nhave not been discussed in the main text. Could this possibly be related to proton delocalization\nasymmetry or subtle interactions within the H-bond network? I also strongly suggest to mark in\nthese figures as well as Fig. 4d which side corresponds to the donor and which to the acceptor,\nas this would help in correctly interpreting the data without looking into them main text.\nReply: We thank the reviewer for the careful observation and helpful suggestion. The small\nfeatures resembling double maxima/minima on the acceptor side in Fig. 4b,c are primarily\nattributed to computational noise due to limited sampling and statistical fluctuations. We improved\nout statistical method to smooth the probability distribution and free-energy profiles. The updated\nresults are shown in Fig. 4a,b.\nAction: We have clearly labeled “donor” and “accepter” in the revised Fig. 4a. The improved\nstatistical method is provided in the SI (page 6).\nComment 7. In Fig. 3e,f, I find it somewhat challenging to clearly distinguish the single and\ndouble lobes that the authors refer to in the discussion. While Fig. 3i,j certainly help clarify the\ninterpretation by quantifying the force differences, perhaps there is a more effective way to visually\nenhance these features in the AFM images.\n−8.5°-L hexamer (Fig. 3c,f) with identical color code. Via comparison, the double-lobe contrast is\nclearly visible for L3-L4, L4-L5 and faintly for L2-L3. We agree that the M-shaped the F (d)\nmin\ncurve provides the more straightforward, quantitative distinction between the NH and N sites (see\nalso Reply (1,2) to Comment 7 from Reviewer #1).\nComment 8. The following sentence is missing a verb: “It is important to note that although the\nexperimentally determined apparent distance between intra- or intermolecular N and NH sites\nsystematically larger or smaller than the DFT-calculated values…”\nReply and Action: We thank the referee for pointing this out. The sentence has been revised to\n“…are systematically …” (page 5).\nComment 9. Is there a reason the authors use \"Supplementary Fig. 3\" instead of a shorter format\nlike \"S3\"?\nReply: We follow the routine of Nature Communication style in the main text.\nComment 10. In the Supplementary Information, the authors use hcp in italics, but in the main\ntext, the formatting is inconsistent. It would be best to standardize the notation throughout the\npaper.\nReply and Action: We have ensured that hcp and fcc sites are consistently denoted using italics\nin the main text (page 5) and SI (page 14, 31).\nComment 11. It may just be my impression, but in Fig. 2e, it looks like the N site and NH site are\nnot exactly chosen at the maxima, but rather slightly next to them. Is this just an optical illusion,\nor could the authors clarify how these positions are precisely defined? Additionally, I noticed a\ndifference in the donor-acceptor distance when comparing Fig. 2e and Fig. S2c. in Fig. 2e, the\nvalue is 460 pm, while in Fig. S2c, the value is 490 pm. Is this difference correct, and if so, what is the\nreason behind it?\nReply: In Fig. 2e we did not select the two extrema of the F (d) curves. AFM ∆f signals have\nmin\nhigher noise than STM z signals, and converting ∆f to force AFM via the Sader-Jarvis method\nfurther amplifies that noise. Because the integration time per point must be limited to minimize\nthe drift during spectroscopy, small differences between neighboring extrema it can be buried in\nnoise, making the true maximum ambiguous.\nIn the revision, we re-measured several single-molecule force spectra with improved signal-\nto-noise ratio, and the overall trend remains consistent. The Δf (d) curves show two clear maxima\nmin\n(points 2 and 17 from left to right), while in the coverted F (d) curves the maximum at the N site\nmin\nis less pronounced. Accordingly, the positions reported in the main text are defined from the\nΔf (d) data, giving a donor–acceptor distance of 465 pm.\nmin\nnoise (occasionally leading to 1–2 points discrepancies), this does not affect the main conclusion\nthat F (N) > F (NH).\nmin min\nAction: We have re-measured the single-molecule force spectroscopy data with improved signal-\nto-noise ratio, and these new data have replaced the previous measurements in both the main text\n(Fig. 2e) and the SI (Supplementary Fig. 4).\nPoint-by-point response to the Reviewer #4’s comments\n“…In my opinion, the experimental data are of high quality, and the observed phenomena are\nboth novel and fascinating. The theoretical calculations, particularly the path-integral\nmolecular dynamics (PIMD) simulations, provide robust support for the experimental findings,\noffering a detailed mechanistic understanding of proton delocalization and quantum effects in\nhydrogen bonds. Overall, this is an excellent piece of work, and I agree with its publication in\nNature Communications. However, several issues need to be addressed to optimize the\nmanuscript further”\nComment 1. Regarding the long-standing debate over the origin of the clear resolution in CO-tip-\nbased chemical bond detection (i.e., whether the bright intermolecular contrast represents a\nhydrogen bond or an artifact caused by CO-tip flexibility), the manuscript mentions that \"line\nspectroscopy measurements can unambiguously determine hydrogen bond directionality,\nproviding a clear resolution to this debate.\" Could the author comment more on this conclusion?\nWhy can local force spectroscopy be more convincing in determining the position of the H atom\nthan the features in constant height imaging like Fig. 2f?\nReply: We address the comments as follows.\n(1) A detailed explanation of how line-FS determines H-bond directionality is provided in Reply\nto Common question 1.\n(2) We consider local FS more reliable than constant-height (CH) AFM imaging for assigning\nproton order for the following reasons:\n(i) AFM imaging is tip-height sensitive. Since Δf(z) is non-monotonic, CH images can\nundergo contrast reversals with small height changes, causing misassignment. FS first\nlocates the Δf minimum and thus sets the correct imaging height.\n(ii) Height corrugation and uneven adsorption can obscure key information. For BI chains on\nAu(111), both BR-AFM and proton-order-resolved AFM images are ambiguous, whereas\nFS clearly reveals the proton order.\n(iii) FS yields observables, e.g., donor–acceptor apparent distances and site-specific\ninteraction strengths, which are not readily available from imaging alone.\nand expect it to work for homonuclear O–H···O, and even for heteronuclear N–H···O bonds\nwith simulation support. AFM imaging remains valuable for speed. A practical strategy is to\nuse FS to determine proton ordering on representative structures, then image at optimized\nheight to survey many cases. In our study, FS established the hexamer ground-state, and AFM\nimaging tracked switching during manipulation (Fig. 5).\nComment 2. This work reports a remarkable observation of symmetry-broken quantum proton\nordering in the ground states of chiral cyclic H-bonds within BI hexamers on Ag(111), despite the\nsymmetric adsorption registry as shown in Fig. 3c-d. To clarify this intriguing result, could the\nauthors comment more on this phenomenon? Could it be a local distortion of this hexamer\nduring the deposition of the molecules or a local imperfection of the Ag surface?\nReply: We address raised questions as follows.\nAFM simulations: During the revision, we simulated AFM images of BI hexamers on Ag(111)\nusing PPAFM code combined with DFT-calculated Hartree potential. The simulated AFM contrast\nnear F shows a nearly sixfold, flower-like pattern with six bright lobes (Fig. 3f; Supplementary\nmin\nFig. 21). In contrast, experimental AFM images at F do not exhibit this symmetry (Fig. 3c),\nmin\ndue to the coexistence of single- and double-lobe contrast, implying the presence of both classical\nand non-classical H-bonds.\nRelation between local distortion and NQEs: PIMD simulation of a constrained hexamer on\nAg(111) reproduces proton delocalization in shortened H-bonds (Fig. 4d and Supplementary Fig.\n27). Notably, a slight shortening of apparent H-bond lengths at positions of double-lobe feature\nhas been experimentally identified (Supplementary Fig. 25). Importantly, this distortion cannot be\nreproduced in DFT calculations that do not account for NQEs (Supplementary Fig. 26). It emerges\nwhen NQEs is included as revealed in our DFT-based PIMD calculations (Fig. 4c) and in prior\ntheoretical studies27.\nExcluding surface imperfections: Local imperfections of the Ag(111) surface can be ruled out\nbecause (i) the symmetry-breaking pattern is repeatedly observed across many hexamers;\n(ii) it\ncan be reversibly tuned by tip manipulation (Fig. 5); and (iii) it is also observed for BI-\nhexamer/ Au(111) (Supplementary Fig. 36).\nAction: Figures have been updated in the main text incorporating DFT and PPAFM\nsimulations\n(Fig. 3), as well as PIMD calculations with classical force field performed at 5 K (Fig. 4). Line-FS\ncharacterization of BI-hexamers on Au(111) have been added into the SI (Supplementary Fig. 36).\nComment 3. The PIMD simulations in Fig. 4 elegantly demonstrate proton delocalization in a gas-\nphase BI dimer (L3–L4). However, the quantum behavior of the full BI hexamer adsorbed on\nAg(111)—where six H-bonds interact cooperatively—might not be fully captured by this\n23\nsimplified dimer model. Could the authors comment more on whether the dimer-based simulations\nare accurate enough to describe the cooperative transfer of protons in a hexamer?\nReply: We thank the reviewer for this insightful and important question. Our responses are as\nfollows:\n(1) The full BI hexamer adsorbed on Ag(111) system contains 594 atoms. PIMD simulations at\nthe DFT level for such a large system are computationally unaffordable. To balance\ncomputational cost and the accuracy of the described quantum properties, we chose a\nrepresentative dimer model to represent and characterize the experimentally observed NQE-\ninduced intermolecular proton delocalization in the hexamer. We acknowledge that this\nsimplification does not account for the concerted motions of multiple protons, but it clearly\nshows the quantum nature of proton delocalization between BI molecules.\n(2) To further assess the validity of the dimer model, we further performed PIMD simulations of\nthe full hexamer adsorbed on Ag(111) using the reaction force field8 within the LAMMPS\nframework7. To reproduce experimental results, the N···N distances of L3-L4 and L4-L5 in\nthe simulation were adjusted to be slightly shorter than the other N···N spacing, as shown in\nSupplementary Fig 27a. The resulting proton probability profiles (Supplementary Fig. 27d,e)\nare in qualitative agreement with the dimer PIMD results in Fig. 4a, which verifies the validity\nof the dimer model, i.e. the simplified dimer model can accurately reflect the fundamental\nquantum behavior of proton delocalization.\nAction: We revised Fig. 4 in the main text to include PIMD simulations of a constrained hexamer\non Ag(111) using the reaction force field (see also Supplementary Fig. 27).\nComment 4. The experimental observations of symmetry-broken H-bonds and their manipulation\n(e.g., Fig. 3, 5) were conducted at 4.8 K, while the PIMD simulations (Fig. 4) were performed at\ntemperature (50 K) impact the agreement between the predicted proton delocalization behavior\nand the experimental results obtained at ultralow temperatures (4.8 K)?\nReply: We provide a detailed response to this question in the Reply to Comment question 3.\nComment 5. The manuscript highlights the critical role of adsorption geometry in stabilizing non-\nclassical H-bonds (e.g., symmetric H-bonds in hexamers vs. classical H-bonds in chains;\nSupplementary Figs. 12–13). Could tailored modifications of adsorption sites (e.g., substrate step\nedges, strain engineering) enable the emergence of symmetric H-bonds even in chain\nconfigurations? And why is the delocalized proton only found in the hexamer?\nReply: We thank the reviewer for this suggestion. We also checked if non-classical H-bonds can\nalso be induced in the chains. To this end, we use tip manipulation to bend a linear chain to\nintroduce possible change of H-bond lengths and primary evidence support our interpretation.\non Ag(110). Nano Lett. 18, 2937-2942 (2018).\nPhys. Rev. B 90, 085421 (2014).\nNano Lett. 16, 1974-1980 (2016).\nby atomic force microscopy. Nat. Commun. 9, 122 (2018).\ntetramer on NaCl(001). J. Chem. Phys. 148, 102329 (2018).\nPhys. 11, 235-239 (2015).\nmodeling at the atomic, meso, and continuum scales. Comput. Phys. Commun. 271, 108171\n(2022).\nHydrocarbon/Water Weak Interactions in the Condensed Phase. J. Phys. Chem. B 122, 4083-\n4092 (2018).\non metal surfaces. Science 377, 315-319 (2022).\nStudy of the H5O2+·Ar Vibrational Predissociation Spectra: Exploring the Putative Role of\nFermi Resonances in the Bridging Proton Fundamentals. J. Phys. Chem. B 112, 321-327\n(2008).\naspartate aminotransferase to facilitate the second half-reaction. Chemical Science 13, 10057-\n10065 (2022).\nChemPhysChem 6, 1738-1741 (2005).\nconducting polymers and liquids. Electrochimica Acta 43, 1281-1288 (1998).\n187-274 (Academic Press, 1987).\nNeutralization–Reionization and ab initio Study. Journal of Mass Spectrometry 31, 1173-1184\n(1996).\nmodulation force spectroscopy. Appl. Phys. Lett. 84, 1801-1803 (2004).\nPorphycene Molecules on a Cu(110) Surface. Phys. Rev. Lett. 111, 246101 (2013).\n25\nbenzimidazoles. Nat. Commun. 3, 1308 (2012).\nScience 264, 1887-1890 (1994).\nserine proteases. Science 264, 1927-1930 (1994).\nInt. Ed. 57, 3888-3908 (2018).\nInstrum 90, 011101 (2019).\n— imidazole on grafoil. Surface Science 191, 547-578 (1987).\nimidazole from inelastic neutron scattering. J. Chem. Phys. 103, 8444-8455 (1995).\nFrom Gate-Opening Phenomenon to Shear-Driven Structural Destabilization. Phys. Rev. Lett.\n113, 215502 (2014).\nZeolitic Imidazolate Framework-71. The Journal of Physical Chemistry Letters 13, 2838-2844\n(2022).\nhydroxyl overlayers on metal surfaces. Phys. Rev. Lett. 104, 066102 (2010).\nWe are grateful to all reviewers for their valuable time, thoughtful suggestions, and for helping\nsubstantially elevate the quality of our work. In this revision, we address the remaining questions\nby providing additional experimental data and simulations, along with further analysis. Figure and\nreference numbering refer to those in the revised manuscript and supplementary information (SI).\nReviewer #1\n“The authors provided extensive new data and explanations and clarified most points.\nI think the addition of the bond-order related contrast of C=N and C-N bonds improved the MS.\nThe DFT calculations showing a flat-lying geometry of the hexamers (Supplementary Figs. 10\nand 17), is an important improvement.”\n“The most important point in my opinion is the assignment of the switching to proton tunneling.\nIt is great that the authors performed experiments with deuterated molecules. Showing a\ndifferent voltage threshold for H and D tunneling would be convincing proof of the switching\nbeing related to proton/deuterium tunneling. This the authors attempt to show in the new Fig. S34.\nHowever, with the few individual reported voltages for switching and the way they are presented\nit is not convincing. This could be improved. It is not clear if that is a proof or not. I think without\ntaking additional data this question can and should be answered before publication. It would be\nnice and a very convincing proof of proton/deuterium tunneling if different voltages thresholds\ncould indeed be extracted from the experiment. But, I would also not object to publication, if\nnot.”\nComment 1 and 2. My points 1-4 refer all to Fig. S34. The voltages of switching should better\n(or at least in addition) be shown against the current at which the switch occurred. Assuming a\nprocess mediated by tunneling electrons, the relation to the tunneling current would be more direct\nthan to the tip heigh. It would also remove possible effects that might be related to the set point\nbeing on a different surface (Ag) and at a different voltage compared to the voltage at which the\nreaction occured with the tip on the molecule, and most important\neffects due to different tips being used. Please plot switching Voltages against the current at which the switch occurred.\nAs a side note: The different used tip heights z, z’, Z make it extremely hard to understand and add\nconfusion. The tip height, z’ seems given wrt. z with seems given wrt. the STM setpoint – one has\nto find different passages of the MS and SI to find that out. Is this here z, or z’ or Z? And what are\nthe differences?\nAssuming an electron mediated process (which the authors do, and I agree with that),\nlikely with\na small yield, the voltage of switching will be stochastic. That is, the individual switching voltage\nwill vary in different individual experiments (therefore I would not name that individual switching\nvoltage “voltage threshold”), even under identical experimental conditions. Moreover, the\nswitching voltage will even depend on the ramp velocity of the voltage. That means:\n2a) The (statistical error of the points in Suppl. Fig. 34e) should be estimated. Ideally this is done\nby repeating switching with identical parameters several times (at least for one condition). I assume\n1\nthat authors have such data. Please provide error bars for the points in Supplementary Fig. 34e\n(and plot against current).\n2b) The parameters of the bias sweeps (including the time for the sweep) should be similar or at\nbest identical for the experiments on H and D molecules. And the parameters should be provided.\nI assume the trend in Figure 34e will be the same as now when plotting threshold Voltages as a\nfunction of current. Which would be more convincing. However, not sure how large the error bars\nwill be. This cannot be inferred with the information I have. I assume the error bars will be of the\nsame order of magnitude as the observed differences in switching voltages between H and D\nmolecules.\nReply: We thank the reviewer for these valuable suggestions and we provide our responses as\nfollows:\n(1) Switching voltage vs. current:\n(i) We display switching voltage versus tunneling current with expanded data points in\nSupplementary Fig. 37e. For both H-BI and D-BI hexamers, the switching voltage\ndecreases as the current increases.\n(ii) Note that in the voltage–current plot, each collective proton-transfer event yields a (V, I)\npair, and different events yield different pairs. Thus, error bars at identical voltage or\ncurrent are hard to define. We therefore plot all individual events for H-BI and D-BI\nhexamers to enable a direct comparison.\n(iii) Clearly, within the same current regime, the switching voltages of D-BI hexamers lie\nsystematically below those of H-BI hexamers.\n(2) Side note on z, z’, Z: This point is addressed in our Reply to Comment 5.\n(3) Statistical error for switching voltage vs. tip height:\n(i) We now provide error bars (mean ± s.d.) for the data points in Supplementary Fig. 37d. For\neach tip height, we repeated the switching manipulation multiple times on the same\nhexamer configurations, +8.5°-L (state 2) or −8.5°-R (state 3), using an identical bias ramp\nvelocity.\n(ii) With error bars included, the overall trend remains unchanged: for both H-BI and D-BI,\nthe switching voltage signifying collective proton/deuteron transfer within the hexamer\ndecreases as the tip height decreases.\n(iii)At the same tip height, the switching voltages of D-BI hexamers are lower than those of\nH-BI hexamers.\n(4) Voltage ramp rate: The bias sweep rate was kept constant at 2 mV/s for all manipulation events.\nThis is now stated explicitly in the figure caption.\n2\n(5) Discussion: Plotting the switching voltage against tunneling current reveals a clear trend: the\nswitching voltage decreases with increasing tunneling current, consistent with an electron-\nmediated process. However, the number of switching events is limited. In our experiments, a\nvery strong tip–molecule interaction (~20 nA in current) must first be established to modify\nthe adsorption geometry; only then can collective proton transfer be triggered. This demanding\nmanipulation protocol restricts the number of events that can be recorded with the same tip.\nNevertheless, in Supplementary Fig. 37d,e, at the same tip height or within same current\nregime, D-BI hexamers switch at lower bias voltages than H-BI hexamers, consistent with the\nlower N-D stretching-mode energy (Ref. 47).\nAction:\n(i). We provided the statistical error bars in Supplementary Fig. 37d.\n(ii). We added the voltage ramp rate into the figure caption of Supplementary Fig. 37.\n(iii). We plotted switching voltage against tunneling current in Supplementary Fig. 37e.\nComment 3. Minor points in Suppl. Fig. 34. Should not the labels “NH” better read “ND”, as\nthese experiments were on deuterated molecules?\nReply and action: We have changed the labels from “NH” to “ND” in the figure to reflect that\nthese experiments were performed on deuterated molecules (Supplementary Fig. 37a-c).\nComment 4. On the y-axis in panel e, the values for 280 and 320 mV seem swapped.\nReply and action: We thank the reviewer for careful reading. The 280 mV and 320 mV labels on\nthe y-axis in panel (e) were indeed interchanged. We have corrected the labels and updated the\nfigure in the revised manuscript (Supplementary Fig. 37d). We apologize for the oversight.\nComment 5. It seems at least three different tip heights, z, z’ and Z in Suppl. Fig. 15 in some other\ninstances (Suppl. Fig. 34) are used. The use of the different z, and z’ (and Z) seems extremely\ninconsistent and unclear to me. I understood from the text z’ is used for spectroscopy, but it is then\nmentioned as tip height for many AFM images in the SI. E.g. in S. Fig. 9, z is used in the panels,\nbut z’ in the caption. Both seem to denote the same thing and have the same values. But in the in the\nmethod section\nit says that z and z’ differ by 30 pm. Please check all instances and make clear what is the difference between z and z’.\nDo you really need both of them, why not subtract the 30 pm and only use z, in all instances?\nReply: We agree tip heights should be reported consistently. Following the reviewer’s suggestion,\nwe have revised all relevant figures in the main text and SI so that AFM images and force spectra\nuse a common reference point, defined with respect to the STM setpoint. We now use only two\nsymbols: z for experimental tip heights (AFM image and force spectroscopy), and z’ for the\nsimulation heights.\nThe former “Z” notation has been removed; (iii) The explicit definition has been added into SI\n(page 5).\nComment 6. I focused on Fig. S34 and reporting the things to be improved only in this Figure\ntook me a while. I did not do that so carefully on all figures. Please check yourself. One other thing\nI did spot: Fig. S7: Delta f scales in c and d are very different compared to the other panels. Please\ndouble check, if correct. If correct, maybe comment on it. Was a different tip used?\nReply: In Fig. S7c,d, the Δf scales indeed differ from the other panels. This most likely stems from\na slightly blunter tip, which increases van der Waals interaction between the tip base and the\nmolecule, leading to larger |Δf|. Nevertheless, this does not significantly affect the imaging quality.\nAction: To avoid confusion, we added a note to the caption of Supplementary Fig. 8, stating that\nimages in panels (c, d) were recorded with a different tip.\nComment 7. I now realized that the method to distinguish N and NH sites is very similar to what\nwas shown in F. Mohn, L. Gross, N. Moll, and G. Meyer, “Imaging the charge distribution within\na single molecule,” Nat. Nano. 7, 227–231, 2012, see Supplementary Figure S2. Compare e.g. with\nFigure S15e. This seems to support the analysis of the current MS. The work should be cited and\nmentioned in this context. I think it also strengthens the assignment. Note that in the cited work,\nthe positions of N and NH groups could be inferred independently by orbital density maps\nmeasured with STM. Moreover, the relation of the AFM contrast, to the KPFM contrast (mainly\nrelated to electrostatic interaction) was measured and discussed.\nReply: We thank the reviewer for pointing us to this reference. In that study, the NH and N sites\nin the core of a naphthalocyanine molecule on NaCl/Cu(111) are clearly distinguished by force\nspectroscopy (Fig. S2 in Ref. 31). Notably, at the tip-sample distance corresponding to the\nminimum of the ∆f(z) spectra, the AFM contrast exhibits twofold symmetry: the inner N sites\nappear brighter than the NH sites (inset of Fig. S2 in Ref. 31). The contrast was attributed to short-\nrange electric fields arising from submolecular charge variations, corroborated by KPFM mapping.\nThese observations are consistent with our results and support our assignment of NH and N sites\nin the imidazole ring using the line-FS approach. Importantly, our work further shows that even\nwhen molecular states are absent near the Fermi level–prohibiting site discrimination by orbital\nimaging–CO-AFM/FS can still discern proton ordering.\nAction: We have cited and discussed the suggested reference (Ref. 31) in the main text to further\nsupport our interpretation (page 5).\n“The authors have substantially revised the manuscript, adding height-dependent AFM/FS data,\nisotope substitution, extended simulations (DFT + PPAFM, PIMD, classical MD), and additional\ncontrols. These additions significantly improve the robustness of the study.\nThe central claim—that AFM line frequency-shift profiles can reveal proton delocalization driven\nby nuclear quantum effects—remains controversial, as electrostatics, adsorption registry, and fast\nhopping are not entirely excluded. Nevertheless, the expanded dataset now provides a consistent\nand technically rigorous case that certain AFM signatures are consistent with non-classical\nhydrogen bonds, and that these signatures cannot be captured by classical models or simple\ngeometric variations. With careful framing, this contribution will be of broad interdisciplinary\ninterest and is suitable for publication in Nature Communications.\nI therefore recommend acceptance after minor revision, contingent on the following\nclarifications and additions”\nComment 1.\nThe manuscript establishes a correlation between AFM/FS signatures and the presence of non-\nclassical H-bonds. However, stating that AFM directly resolves or discerns proton delocalization\nseems too strong. I encourage the authors to moderate the wording, using terms such as ‘consistent\nwith’ or ‘compatible with’ rather than categorical claims.\nSections where the language should be softened\nAbstract (p.2): “…BR-AFM/AFS to directly probe the proton ordering…”\nSuggest: “…BR-AFM/AFS to probe signatures of proton ordering…”\nAbstract (p.2): “…we determine H-bonding directionality and detect signature of quantum proton\ndelocalization.”\nSuggest: “…we extract information consistent with H-bonding directionality and signatures\ncompatible with quantum proton delocalization.”\nResults, “Elucidating hydrogen bond directionality…” (p.4): “…FS distinguishes N from NH\nsites…”\nSuggest: “…FS provides contrast consistent with differences between N and NH sites…”\nResults, hexamer analysis (p.6):\n“…Fmin(d) curves reveal a clockwise N−H···N bond directionality…”\nSuggest: “…Fmin(d) curves show contrast compatible with a clockwise N−H···N bonding\nscheme…”\nSummary (p.10):\n“…a line-FS based method that resolves H-bonding directionality…”\nSuggest: “…a line-FS based method that provides evidence consistent with H-bonding\ndirectionality…”\nmanuscript to avoid overstatement. Corresponding changes are highlighted on pages 2,4,6 and 11.\nComment 2. Some concerns that still need clarification from the authors. Quantitative robustness\nof the classification criterion:\nReport the number of spectra averaged per bond and provide error bars for |ΔF_min|.\nSummarize the fraction of bonds classified as non-classical across all measured.\nReply: We address reviewer’s concerns as follows:\n(i) Number of spectra averaged per bond: For each line-FS across a given H-bond, we record a\nsingle f(z) spectrum at each point along the path, using a sufficiently long integration time (~1 min\nper spectrum) to ensure a high signal-to-noise ratio. The dominant source of uncertainty is thermal\ndrift, which can slightly offset the probed sites relative to the N and NH positions. To minimize\ndrift, once a hexamer is located, we keep the tip tunneling over the same region for ≥12 h. Under\nthese conditions, no significant drift is observed after each line-FS measurement.\n(ii) Error bars for |∆F |: Systematically repeating every line-FS at each individual H-bond (e.g.,\nmin\nthree repeats per bond) would require ~9 h of continuous measurement per hexamer, which is\nexperimentally very demanding and practically difficult to achieve. Therefore, not every H-bond\nin every hexamer was repeatedly measured. Nevertheless, in selected cases, we did perform\nrepeated line-FS measurements on the same hexamer under the same tip condition to estimate\nuncertainties. The resulting error bars are small and do not affect our conclusions, thereby\nvalidating the robustness of our H-bond classification.\n(iii) Fraction of non-classical bonds across measured hexamers: We went through the data of all\nhexamers characterized by line-FS or proton-order sensitive imaging (56 cases). In every case, the\nnumber of non-classical H-bonds per hexamer was either 1, 2 or 3 out of 6 and the mean fraction\nof non-classical H-bonds is approximately 40%.\nAction: (i) We now include a section in the SI detailing the line-FS measurements (see “Line\nforce spectroscopy measurements on H-bonds” on page 5); (ii) For clarity, we now include\nrepresentative datasets with error bars in the SI (see Supplementary Fig. 17c,d), demonstrating the\nhigh stability and reproducibility of our measurements; (iii) The mean fraction of non-classical H-\nbonds in all measured hexamers has been added into the main text (page 7).\nComment 3. Height-dependent analysis: Present Δf(z) traces at N vs NH with residuals from your\nfits, to show explicitly the z-range over which the classification is stable. Indicate the uncertainty\nfrom the Sader–Jarvis conversion.\nReply: Following reviewer’s suggestion, we provided force-spectroscopy data processing details\nin the SI (Supplementary Fig. 17a,b), presenting the raw Δ𝑓(𝑧) data, the smoothed curves, and the\ncorresponding residuals. As shown in Supplementary Fig. 17a, the Δ𝑓 measurements exhibit a\nvery high signal-to-noise ratio, with residuals around the minima < 0.02 Hz.\nin the input Δ𝑓 data, particularly at short tip-sample distances where Δ𝑓(𝑧) changes quickly. To\nalleviate the resulting noise in the 𝐹(𝑧) curves, we applied a mild smoothing procedure and show\nboth the unsmoothed and smoothed data together with their residuals. As shown in Supplementary\nFig 17b, although the residuals reach ~3 pN at the closest approach, they remain ~ 1 pN in the\nvicinity of the curve minima–the region relevant to our analysis. This level of uncertainty does not\naffect the determination of the H-bond directionality or the classification into quantum-delocalized\nversus classical H-bonds.\nAction: We included the related analysis and discussion into SI (Supplementary Fig. 17)\nComment 4. Internal charge effects: Internal charge redistribution within a single BI molecule is\nexpected to distort AFM contrast patterns directly impacting into bond-length analysis from\nexperimental images. Did the authors explicitly check for such distortions and compare with bond-\nlength analysis obtained from calculations? Please clarify.\nReply: We agree that internal charge redistribution within a single BI molecule can influence\nbond-resolved AFM images via electrostatic interactions with a CO-tip. To assess possible contrast\ndistortions from these effects, we performed additional AFM simulations for a hexamer on Ag(111)\nby systematically varying the tip charge. The main findings are as follows.\nAs shown in Supplementary Fig. 14, a negatively charged CO-tip clearly resolves the five-\nmembered (imidazole) and six-membered (benzene) rings, consistent with experiment. In contrast,\nneutral or positively charged tips yield distorted, contracted five-membered rings, deviating from\nthe experimental observation. Bond-length analysis across tip-charge states shows that, although\nabsolute apparent bond lengths shift with tip charge, the relative trend—C–N bonds being longer\nthan C=N bonds—remains robust, being coherent with a classical chiral cyclic H-bonding scheme.\nOverall, these results indicate that electrostatics are essential for reliable AFM simulations. In our\ncase, assigning a negative charge to the CO-tip reproduces both proton-sensitive contrast and\nbond-resolved AFM images.\nAction: (i) We added a detailed discussion of tip-charge effects on bond-length analysis of\nsimulated AFM images in the SI (Supplementary Fig. 14); (ii) we cited J. van der Lit et al., Phys.\nRev. Lett. 116, 096102, (2016), which specifically addresses the role of electrostatics in high-\nresolution AFM imaging in the main text (Ref. 38); (iii) a minor amendment has been added into\nthe main text on page 6.\nComments: The authors have satisfactorily addressed the points I raised in my previous report.\nAs I am not a specialist in the specific experimental methodology, I am not in a position to\ncomment in detail on some of the technical issues discussed in response to other reviewers’\ncomments and will leave their evaluation to experts in the field. From my side, I have no further\ncomments.\nReply: We are pleased that our revisions have satisfactorily addressed Reviewer 3’s concerns, and\nwe appreciate the reviewer’s expert evaluation.\nReviewer #4\nComments: The authors have thoroughly addressed all comments raised in this round of review.\nThey provided further/newer data where necessary and added relevant explanations where required.\nThe paper reads more clearly, and the analysis of the results is more convincing. This study\nadvances our understanding of the properties of cyclic hydrogen bonds. I therefore recommend\npublication of “Anomalous symmetry breaking in proton ordering of surface-confined cyclic\nhydrogen bonds” in Nature Communications.\nReply: We are glad that the added data and explanations made the manuscript clearer and the\nanalysis more convincing. We appreciate Reviewer 4’s positive evaluation and recommendation\nfor publication.\n8"}]

Decision letter