Elemental Identification by Combining Atomic Force Microscopy and Kelvin Probe Force Microscopy

There are currently no experimental techniques that combine atomic-resolution imaging with elemental sensitivity and chemical fingerprinting on single molecules. The advent of using molecular-modified tips in noncontact atomic force microscopy (nc-AFM) has made it possible to image (planar) molecules with atomic resolution. However, the mechanisms responsible for elemental contrast with passivated tips are not fully understood. Here, we investigate elemental contrast by carrying out both nc-AFM and Kelvin probe force microscopy (KPFM) experiments on epitaxial monolayer hexagonal boron nitride (hBN) on Ir(111). The hBN overlayer is inert, and the in-plane bonds connecting nearest-neighbor boron and nitrogen atoms possess strong covalent character and a bond length of only ∼1.45 Å. Nevertheless, constant-height maps of both the frequency shift Δf and the local contact potential difference exhibit striking sublattice asymmetry. We match the different atomic sites with the observed contrast by comparison with nc-AFM image simulations based on the density functional theory optimized hBN/Ir(111) geometry, which yields detailed information on the origin of the atomic-scale contrast.

A tomic-resolution microscopies are key enabling techniques in modern materials research. These techniques include scanning probe microscopies (scanning tunneling microscopy, STM, and noncontact atomic force microscopy, nc-AFM) 1−3 as well as high-resolution transmission electron microscopy (e.g., scanning transmission electron microscopy, STEM). 4 Identifying different elements is more challenging, especially based purely on experimental results. 5−9 Modern, aberration-corrected electron microscopes make it possible to chemically fingerprint different elements through atomic-resolution electron energy loss spectroscopy (EELS) 8−11 or Z-contrast in the annular dark field imaging mode. 12 These techniques, when coupled with ab initio calculations, can also yield information on chemical bonding configuration down to the single-atom level. 10,11 However, electron microscopy has demonstrated chemical sensitivity only on solid state materials. Due to the high-energy electron beam, STEM has not yet reached atomic structural or chemical resolution on more sensitive structures such as small organic molecules. On the other hand, scanning probe microscopy can be used to image small molecules with atomic resolution. 13,14 When operated in the frequency modulation mode, 15 nc-AFM measures atomic-scale forces between the tip on an oscillating cantilever and the sample surface through changes of the resonance frequency (Δf) of the cantilever. nc-AFM can yield atomic resolution, which has been amply demonstrated on elemental semiconductors as well as on heteroatomic surfaces of compound semiconductors and polar insulators such as alkali halides and oxides. 1,3,6,16−18 On such heteroatomic surfaces, typically only one type of atom is actually imaged with a given tip, because the polar nature of the compounds results in a strong variation in the short-range forces above the negatively and positively charged atoms. 5, 19 Additional information can be gained from force spectroscopy, i.e., the measurement of Δf as a function of tip−sample distance z. In a seminal contribution, Sugimoto et al. identified different atomic species in a disordered surface alloy on Si(111) by comparing the maximum attractive forces on different lattice sites. 7 This methodology has recently been extended to estimating electronegativities of surface atoms. 20 In addition to the force channel, another avenue for achieving elemental contrast is to measure the local contact potential difference (LCPD) by Kelvin probe force microscopy (KPFM). 21−27 Nevertheless, contrast between different elements on the surface has thus far only been observed for highly polar or reactive surfaces.
Atomic-resolution nc-AFM studies can be extended to molecular systems through chemical passivation of the tip apex, e.g., by controlled pick-up of a single carbon monoxide (CO) molecule. 13,14,28−36 With these tips, it is possible to enter a regime where the tip−sample interaction is dominated by the Pauli repulsion between the last atom of the tip and the sample atom directly under it. 28,37−39 In addition to molecules, this technique has been used to measure atomic positions and surface corrugations of two-dimensional materials (e.g., graphene and hexagonal boron nitride). 32,40−44 Achieving chemical sensitivity with passivated tips is a more delicate issue, as bond formation with the sample similar to reactive tips is suppressed. Still, changes in the total electron density and electrostatic forces on different atoms are expected to contribute to the image contrast. 18,37,45−48 Consequently, passivated tips enabled elemental contrast also on molecular systems, in both Δf 29,49−52 and LCPD. 50,53 However, cross-talk between perceived elemental contrast and nonplanar topography, edge effects, and resulting image distortions due to flexibility of the tip apex make systematic studies in molecular systems challenging. 34,35,51,52,54−56 At present, chemical fingerprinting by scanning probe microscopy remains an elusive goal.
Here, we demonstrate elemental contrast in both Δf and LCPD on a covalently bonded system, monolayer hexagonal boron nitride (hBN). This system is well-defined on the atomic level and does not suffer from topographic corrugation or edge effects. We employ nc-AFM with CO-functionalized tips 28 to investigate the atomic-scale contrast on epitaxial hBN on Ir(111). 57 Despite the mostly covalent character of the B−N bond 58,59 and a nearest-neighbor distance of only ∼1.45 Å, constant-height maps of both Δf and LCPD acquired over hBN/Ir(111) exhibit striking sublattice asymmetry. nc-AFM image simulations based on the density functional theory (DFT)-optimized hBN/Ir(111) geometry allow us to match the two distinct atomic sites with the boron and the nitrogen sublattices. Studies on such clean model systems are essential to shine light on the origin of atomic-scale contrast in nc-AFM on surfaces and molecular systems. Figure 1a shows the DFT-optimized structure of our model surface, monolayer hBN on Ir(111) (see Methods for computational details). The lattice mismatch between hBN and the iridium surface results in a periodic variation of the stacking between boron and nitrogen and substrate atoms, thus creating a moirésuperstructure with a periodicity of approximately 30 Å. The different registries dictate the interaction strength between the hBN layer and the metallic surface, giving rise to a structural corrugation of the overlayer of ∼1.3 Å, as well as a modulation of the local work function within the moiréunit cell. 57 However, the regions surrounding the moirédepressions are essentially flat. Figure 1b is an atomically resolved STM image of hBN grown on Ir(111) by chemical vapor deposition (see Methods for experimental details), acquired with a CO-passivated tip prepared on a Cu(111) surface (CO/Cu tip) (all data shown here were acquired with CO/Cu tips). The topography highlights the depressions of the superstructure and alignment of the hBN lattice with the moiréunit cell, in good agreement with the DFT calculations and previous studies. 57,60 Two constant-height nc-AFM images, recorded at different tip− Similarly, the atomic contrast comprises a hexagonal lattice of top sites with less negative Δf. Both contrasts seem to reverse at small tip−sample distances (panel d, Δz = 0.7 Å), where the moirédepressions exhibit more negative Δf, and at the atomic scale, a honeycomb lattice of less negative Δf appears.

RESULTS AND DISCUSSION
In the limit of small amplitudes, for which the frequency shift is given by Δf = −f 0 /2k. ∂ z F z (F z is the vertical component of the tip−sample force 1 and f 0 and k are the resonance frequency and stiffness of the cantilever, respectively), more negative Δf values can be interpreted as more attractive tip−sample forces. Then, the moirécontrast reversal is easily explained: At large tip−sample distance, the depressions are less attractive than the surrounding regions, while at close distances they are less repulsive [see Supporting Information (SI) for comparison of the frequency shift as a function of tip−sample distance between the moirédepression and the surrounding region].
Understanding the atomic-scale features requires a more careful inspection, as a mere reversal from attractive to repulsive contrast is not expected for CO-passivated tips. 32 Figure 2a−h show a series of distance-dependent constant-height nc-AFM images, covering a z-range of 1.6 Å. At very large distances (panel a), no atomic contrast is observed, but only the longrange moirécorrugation. Approaching the surface (panel b), the hexagonal pattern of less negative Δf appears, similarly to Figure 1c. On further reducing the tip−sample distance (panels c−e), every second hollow site starts to exhibit less negative Δf as well, such that a lattice with three distinct sites forms. At very close distances (panels f−h), the intensities of the two sites with the less negative Δf approach each other and equalize. Because the two sites have different sizes, this gives rise to a hexagonal lattice of repulsive triangles, which we interpret as the atomic lattice of the hBN. Thus, the sites of less negative Δf observed at large distances ( Figure 2b) do not correspond to the hollow sites of the atomic hBN lattice appearing less attractive but one of the two top sites appearing more repulsive. We have reproduced this evolution of the atomic contrast with several macroscopically different CO tips prepared on Cu(111), as well as with CO tips prepared on Ir(111) [see SI for images acquired with a CO tip prepared on Ir(111)]. The observed contrast is thus intrinsic and not related to tip artifacts.
The above observations reveal a clear asymmetry between the boron and nitrogen sublattices in nc-AFM images of hBN. We can quantify this asymmetry by measuring the frequency shift as a function of tip−sample distance [Δf(z)], as depicted in Figures 2i,j. The graph in Figure 2j shows Δf(z) spectra for the three distinct sites, each of which is the average of three equivalent locations as marked in Figure 2i (see SI for the individual spectra). In Figure 2k, we plot the difference (black) between the two inequivalent top sites (green, red) of the hBN lattice and the corresponding force difference (magenta) recovered from Δf via the Sader−Jarvis method. 61 Taking the adjacent-averaged force data, the difference between the boron and nitrogen sublattice at typical imaging distances is no more than 10 pN. Note that the Δf(z) approach curves go to smaller tip−sample distances than the images in Figure 2a−h, but we observed instabilities of the tip−sample junction when imaging at such small distances.
In order to match the two sublattices with the boron and nitrogen atoms, we simulate nc-AFM images using the MechAFM code, 62 which is based on a molecular mechanics

ACS Nano
Article model taking into account the flexibility of the CO molecule at the tip apex 34,35,41 and the atomic coordinates from the DFToptimized hBN/Ir(111) structure (see Methods for details of the nc-AFM simulations). Recent studies have shown the importance of electrostatic forces in understanding AFM image contrast, 46,47 and thus they need to be taken into account as well. However, there are contradicting reports regarding the charge associated with CO-passivated metal tips. 18,35,45,46,48 One issue is that ab initio simulations of such tips usually focus on unrealistically small tip models, often containing none 28,37,39 or only a few metal atoms 30,31,38,48 in addition to the CO molecule. We address this problem by employing DFT to simulate more realistic tip models, which also contain a metallic surface to account for the bulk tip 56 (see Methods for computational details).
In Figure 3a, we plot the effective long-range electrostatic field of the tip extracted from our DFT calculations for several CO/Cu and CO/Ir tip models. The different models are grouped into four classes: (i) purely metal, pyramidal clusters ("cluster tips"), (ii) metal clusters on a metal surface ("surface tips"), (iii) CO on a metal cluster ("CO-cluster tips"), and (iv) CO on a metal cluster on a metal surface ("CO-surface tips"). An example tip and its actual electric field is shown for each of the four classes in Figure 3b−e (see SI for the structures of all calculated tips). In order to clearly demonstrate the effect of the CO molecule and Cu surface on the electric field, all example tips have the same 10-atom Cu cluster. The cluster tip in Figure  3b exhibits a positive electric field, as is commonly assumed for metallic tips. 18,45,48 However, Figure 3a indicates that within DFT, not just the magnitude but also the sign of this field depend on the precise cluster geometry. For surface tips, we find consistently a positive electric field at the apex, as demonstrated in Figure 3d. Adding a CO to a Cu cluster tip, as shown in Figure 3c, results in a negative electric field at the

ACS Nano
Article apex, in agreement with previous DFT calculations for a COcluster tip. 48 Indeed, we find that most of the CO-cluster tips exhibit a negative effective electric field. In contrast, Figure 3a reveals that within the most realistic subset of CO-surface tips, all but one exhibit a positive effective electric field. This positive electric field is due to the bulk metallic tip and results from the Smoluchowski effect, and this is the dominant long-range contribution. We observe a negative field at around 2 Å from the oxygen apex (see Figure 3e), 48 but the extent of this negative component is exaggerated for cluster tips as seen in Figure 3c, due to the multiple exposed apexes and their unphysical impact on the Smoluchowski effect. In the more realistic CO-surface tip models used here, the field rapidly becomes positive, as shown in Figure 3e. The electric field for the CO-surface tip shown in Figure 3e agrees qualitatively well with the semiempirical CO tip model introduced by Ellner et al., 48 which was fitted to reproduce the nc-AFM contrast of a Cl vacancy in NaCl. However, in their model the negative electric field close to the oxygen extends more into the vacuum compared to our CO-surface tips before it becomes positive. Our finding applies to both Cu and Ir tips, suggesting that this is a general feature of CO/metal tips. Calculation of the dipole associated with the whole tip instead of the effective electric field yields qualitatively identical results, and the conclusions do not depend on which method is used (see SI for a plot of the dipole for all simulated tip models). Hence, the simplest reasonable approximation is to consider the CO tip as carrying a positive partial charge. This charge then interacts with the electrostatic potential of the hBN/Ir(111) sample obtained from DFT, thus accounting for the electrostatic forces in our nc-AFM simulations. Despite its simplicity, such an approach has previously shown to yield excellent agreement with experimental nc-AFM images. 14,47 Three large-scale, simulated nc-AFM images of hBN/Ir(111) at different tip−sample distances and for a positively charged CO tip (q tip = 0.5 e) are shown in Figure 4a. Surprisingly, the image simulations suggest that the boron sublattice appears first in the constant-height images, with the nitrogens becoming visible only as the tip further approaches the hBN layer.
For comparison, we also simulated nc-AFM images with a neutral (q tip = 0.0e) and negatively charged (q tip = −0.5e) tip, as shown in Figure 4b and c, respectively. The main observations are as follows: Even when neglecting electrostatic interactions as for the neutral tip, there is a slight asymmetry between the boron and nitrogen sublattice, with the former appearing first at large tip−sample distances. This atomic contrast is thus due to short-range interactions, which are described by Lennard-Jonestype pair potentials in the nc-AFM simulation, and is in agreement with the larger atomic radius of boron compared with nitrogen, 63 which causes an earlier onset of repulsive forces. This asymmetry then becomes more pronounced with a positively charged tip, which also causes strong distortions of the atomic honeycomb lattice at small tip−sample distances, such that it appears as a hexagonal lattice of repulsive triangles. This contrast is found neither for the neutral nor for the negatively charged tip. On the contrary, the negatively charged tip yields a nearly perfect honeycomb lattice at such small distances. At large tip−sample distances, it shows the nitrogens first and, importantly, also causes a reversal of the moireć ontrast, which is not observed for the other two tips. Thus, by far the best agreement with the experimental data is achieved

ACS Nano
Article for the positively charged tip (see SI for image simulations at additional heights).
The fact that the boron sublattice is observed first is rather counterintuitive: As repulsive atomic contrast is dominated by the Pauli repulsion, one would expect the element exhibiting the higher charge density to appear first upon approaching the surface. Figure 4d shows constant-height slices of the electron density from our DFT calculations, at the heights of the tip's oxygen atom in the image simulations. While the slice at largest distance has significant contributions of the underlying Ir substrate (note also the very small scale of the changes in the electron density), at the two smaller tip−sample distances, the nitrogens clearly exhibit a higher electron density. Thus, our nc-AFM images contradict the conventional interpretation, which highlights the limitations of relying only on the Δf channel. 24 Figure 5a−h display a full series of simulated nc-AFM images for the positively charged CO tip, covering a total range of tip− sample distances of 1.6 Å, the same as in the experimental series (Figure 2a−h). The detailed evolution of the moiréand atomic contrast is in excellent agreement: The theoretical images reproduce all three observed atomic contrasts, as well as the relative range of tip−sample distances at which they appear. In addition, the simulation indicates that only at closest tip− sample distances (Figure 3i) is the contrast formation strongly influenced by the bending of the CO at the tip apex, while this effect is negligible for images taken farther away. It is noteworthy that even though it is well-established that bulk metallic tips carry a positive electric field, 48,64 our finding that it also influences the atomic contrast with CO-terminated tips is surprising with respect to previous studies. 48 In fact, based on the study of the nc-AFM contrast of NaCl and a Cl vacancy, it was shown that the more complicated multipole character of the CO-tip electric field can be relevant in more strongly ionic systems 48 (see SI for details). We have tested this in the case of hBN by carrying out additional image simulations. We used a modified version of the probe−particle model, which calculates the electrostatic forces based on a tip that consists of a positive dipole accounting for the bulk metallic tip and a small negative quadrupole moment for the CO (see SI for details of the nc-AFM simulations). These simulations qualitatively agree with our simpler model and also suggest that the boron lattice is imaged first at larger tip−sample distances (see SI for the simulated images and additional discussion). The fact that already the simple approximation of a positive point charge yields excellent agreement with experiments in the case of hBN is most likely due to the different bonding character compared to NaCl. In the former, the bonds have only small polar character, 59 and thus the electric field associated with the partial atomic charges extends less into the vacuum (in the case of monolayer hBN/Ir(111), the partial atomic charges are further screened by the metal). The electrostatic interaction is then dominated by the long-range positive field of the tip, while its negative part close to the oxygen is negligible. For the ionic lattice of NaCl and the vacancy, the electric field extends much more into the vacuum (for bilayer NaCl there is also less screening compared to monolayer hBN), and thus the interaction with the oxygen's localized negative electric field becomes significant. In addition to these general trends, the effective strengths of the tip dipole and quadrupole moments can vary somewhat in different experiments, which further contributes to the observed contrasts. This suggests that it is important to consider more realistic tip models for some systems and that, in particular, cluster-based models should be used with caution.
To validate the interpretation of boron atoms appearing at larger tip−sample distances than the nitrogens and to provide an additional data channel for elemental identification, we have also performed atomically resolved KPFM experiments. We used a CO-passivated tip to acquire a set of Δf(V) spectra 65 within the area marked with a red square in Figure 6a and extracted the voltage corresponding to the maximum of the Δf(V) parabola, as shown in the example spectra in Figure 6d (see Methods for experimental details). The resulting LCPD map is shown in Figure 6c. In addition, Figure 6b displays a simultaneously recorded nc-AFM image that allows comparison with the experiments shown in Figure 2 and the nc-AFM simulations in Figure 5. The LCPD map shows atomic resolution, with a hexagonal pattern of regions of more negative LCPD values. The simultaneously recorded nc-AFM image allows us to match these regions with the sublattice appearing at larger distances in the nc-AFM images, i.e., the boron atoms according to the image simulations. Importantly, Figure 6b also shows that the elemental contrast in LCPD is achieved at tip−sample distances that correspond to only weak atomic contrast in nc-AFM images, without the necessity to approach the sample into a regime where the signal is influenced by CO bending, tip−sample junction instabilities, or Δf(V) spectra deviating from their expected parabolic shape. 53 At these relatively large tip−sample distances, the LCPD contrast is predominantly governed by the vertical component of the electric field of the sample (E z ). 65 In this approximation, a more positive electric fieldcaused, for example, by partial positive chargesresults in a more negative LCPD. Figure 6e is a map of E z over the entire moiréunit cell, calculated from the Hartree potential of our hBN/Ir(111) DFT simulations. The constant-height slice is taken at 3.8 Å above the mean adsorption height of the hBN layer, i.e., the height of the CO tip's oxygen atom from the nc-AFM image simulation in Figure  5c, which shows good agreement in terms of atomic contrast to Figure 6b. The reduced work function at the depressions of the moiré5 7 results in a long-range modulation of E z , with the depressions exhibiting more positive values. Figure 6f is a zoom-in at the region marked by the red square in Figure 6e, corresponding to the region of the experimental LCPD map (see Figure 6a). The theoretical map shows similar elemental contrast to the experimental one, allowing us to identify the regions of more negative LCPD with the boron sublattice. This

ACS Nano
Article finding is consistent with the interpretation of the nc-AFM images, thus confirming our assignment for the elements on the two sublattices.

CONCLUSIONS
We have demonstrated elemental contrast in nc-AFM and KPFM images on monolayer hexagonal boron nitride. Using a passivated tip and a covalently bonded, inert model surface, our results expand the limits of elemental identification with stateof-the-art atomic force microscopy. Most importantly, our method of combining constant-height Δf images and KPFM maps presents a robust method to identify different chemical species using nc-AFM.

METHODS
Experimental Procedures. Monolayers of hBN on Ir(111) were grown by low-pressure high-temperature chemical vapor deposition of borazine (B 3 N 3 H 6 ), as described in detail in ref 57.
All subsequent nc-AFM and STM measurements were carried out in a Createc LT-STM/AFM, operated at ultrahigh vacuum (UHV) and at a temperature of 5 K. The microscope was equipped with a qPlus 66 tuning fork sensor, which had a resonance frequency f 0 of ∼30.68 kHz, a quality factor Q of ∼98k, and a stiffness k of ∼1.8 kN/ m. For tip functionalization, CO was dosed from a leak valve attached to the UHV system onto the cold surface, either Cu(111) or hBN/ Ir(111). CO was picked up from Cu(111) as described previously. 67 After successful tip passivation, the Cu(111) sample was exchanged for the hBN/Ir(111) sample to carry out the nc-AFM measurements. 32 On Ir(111), CO pickup cannot be carried out as well-controlled as on Cu(111) because of the strong binding of CO to the Ir surface.
Instead, CO pickup was achieved by scanning an Ir(111) surface with high CO coverage at relatively harsh feedback parameters, e.g., ∼10 mV sample bias voltage and several tens of nA tunneling current. This usually led to the transfer of CO from the surface to the tip apex after some scanning time.
All nc-AFM images and KPFM maps were recorded in the constantheight mode with the z-feedback loop disabled and with a tuning fork oscillation amplitude of 50 pm. The sample bias voltage for nc-AFM images was 0 V. KPFM measurements were performed by collecting Δf(V) spectra on a 32 × 32 grid, 65 where each spectrum took 5 s. The spectra were fitted with second-order polynomials to yield the voltage corresponding to the maxima of the Δf(V) parabolas. This voltage minimizes the electrostatic tip−sample forces and is plotted in the LCPD map. Unless stated otherwise, all nc-AFM and KPFM data are unfiltered raw data.
DFT Calculations of hBN/Ir(111). We used the CP2K software package, 68 in particular the QuickStep module, 69 for the DFT calculations of the structure of hBN/Ir(111). The vdW-DF2-B86r approximation 70,71 from the LibXC library 72 was employed for the exchange−correlation term in the Kohn−Sham scheme. The Gaussian plane wave method 73 was used to solve the electronic structure selfconsistently, where the basis set to expand the wave functions was DZVP-MOLOPT-SR-GTH 74 and a cutoff energy of 700 Ry was used for expanding the density in plane waves, with five grids and a relative cutoff of 70 Ry. The Goedecker−Teter−Hutter (GTH) type pseudopotentials 75 were employed, with 17 valence electrons in Ir. Only the Γ point was included in the reciprocal space, without further sampling of the first Brillouin zone. The Fermi−Dirac broadening of the occupation numbers with a "temperature" of 300 K was used. We used a 12 × 12-on-11 × 11 unit cell, in agreement with previous experimental results, 57,60 and four layers of Ir(111). Due to the weak

ACS Nano
Article interactions and strong requirement for precision, the convergency criterion on the maximum force on any ion was set to 0.0514 meV/Å.
DFT Calculations of Different Tip Models. First-principles calculations of the different tip models were performed using the periodic plane-wave basis VASP code 76,77 implementing the spinpolarized DFT. To accurately include van der Waals interactions in this system, we used the optB86b-vdW-DF functional, 78−80 selected based on previous work showing that it provides a sufficiently accurate description for all subsystems involved in the measurement. Projected augmented wave potentials were used to describe the core electrons, 81 with a kinetic energy cutoff of 550 eV (with PREC = accurate). All calculations were performed with dipole correction to account for spurious electrostatic interactions between neighboring cells. Systematic k-point convergence was checked for all periodic systems, with sampling chosen according to system size. This approach converged the total energy of all the systems to the order of meV. The properties of the bulk and surface of Cu, Ir, the isolated structure of CO, and its adsorption on Cu and Ir were carefully checked within this methodology, and excellent agreement was achieved with experiments. A vacuum gap of at least 1.5 nm was used in general, with larger gaps to study the tip potentials. All systems considered were relaxed until the atomic forces were less than 0.01 eV/Å. The tips' electrostatic field was calculated by exploring the electrostatic potential from the DFT simulations. The average electrostatic potential in x−y was calculated as a function of z for each tip model. Beyond the tip apex, this rapidly converges to a constant slope (which reflects the compensation of the electrostatic interactions across the periodic cells and is directly correlated with the tip dipole 82 ). While the absolute value of this slope is not particularly meaningful, it gives a clear relative measure of the effective electric field at long range. For comparison, the conventionally calculated total dipole in the z-direction for each model is shown in the SI, and it shows the same trends across different classes of tip.
nc-AFM Image Simulations. nc-AFM simulations were based on the probe particle model, 34,35,41 widely used to model the nc-AFM imaging process with functionalized tips, 14 as implemented in the MechAFM code. 62 The probe consists of a fixed C atom holder connected to an O atom restrained in the xy-plane by a harmonic spring. The tip−sample interactions were calculated by placing the tip in several locations above the DFT-optimized structure of hBN/ Ir(111) and relaxing the tip O termination. The interatomic interactions between O and the atomic species of the sample are described by Lennard-Jones potentials. For C and O atoms, we used the CHARMM force field parameters. 83 For B and N atoms, the parameters from Hilder et al., 84 derived for the interaction of boron nitride nanotubes with water (H 2 O), were used. Parameters for pair potentials were obtained using arithmetic mixing rules from the atomic ones. In addition, the MechAFM code allows for inclusion of electrostatic forces by assigning a charge to the O atom and letting it interact with the Hartree potential of the sample as obtained from DFT. Frequency shift images were calculated from the tip−sample interaction maps using the method in ref 1, assuming for the O atom a harmonic spring stiffness of 0.5 N/m and a charge of 0.5, 0, or −0.5e. The probe particle model including more complicated tip electrostatics is discussed in the SI.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsnano.7b08997.
Additional nc-AFM data and simulations and structure of all calculated tips (PDF)