Andrea Pérez-Villa

Thomas Georgelin

Jean-Franc ¸ois Lambert

Marie-Christine Maurel

Franc ¸ois Guyot

A Marco Saitta
email: marco.saitta@sorbonne-universite.fr

Fabio Pietrucci
email: fabio.pietrucci@sorbonne-universite.fr

A common precursor to both prebiotic and biological pathways to RNA nucleotides

Understanding the mechanism of spontaneous formation of ribonucleotides under realistic prebiotic conditions is a key open issue of origins-of-life research. In cells, de novo and salvage nucleotide enzymatic synthesis combines 5-phospho-α-D-ribose-1-diphosphate (α-PRPP) and nucleobases. Interestingly, these reactants are also known as prebiotically plausible compounds. Combining ab initio simulations with mass spectrometry experiments, we compellingly demonstrate that nucleobases and α-PRPP spontaneously combine, through the same facile mechanism, forming both purine and pyrimidine ribonucleotides, under mild hydrothermal conditions. Surprisingly, this mechanism is very similar to the biological one, and yields ribonucleotides with the same anomeric carbon chirality as in biological systems. These results suggest that natural selection might have optimizedthrough enzymes -a pre-existing ribonucleotide formation mechanism, carrying it forward to modern life forms.

Introduction

Exploring the synthesis and functionality of nucleic acid molecules is one of the key issues in the origins-of-life research. Since the late 80's, scientists have proposed (1) the possibility of a 'RNA world' as the precursor of current biochemistry, a hypothesis strongly supported by the high versatility of RNA molecules to perform different biochemical processes, such as catalysis, gene regulation, storage of information and self-replication, that are fundamental features involved in the development and propagation of life. A plausible 'RNA World' implies sufficiently successful synthesis of the RNA monomers, namely ribonucleotides (or simply nucleotides), followed by their accumulation, their subsequent polymerization and their further propagation and persistence in the environment of occurrence (2,3). All these steps are affected by a great number of factors (temperature, presence of catalytic surfaces, pH, ionic strength, and so forth), and their plausibility is mainly determined by the thermodynamics and kinetics of the chemical reactions involved.

In the last decades, many efforts have been made in addressing some of these questions.

Nucleotides are composed of three chemical subunits: a ribose in its β-furanose conformation, a nucleobase (purine or pyrimidine) and a phosphate group at the 5' position of the sugar.

One of the most challenging issues has been the synthesis of both purine and pyrimidine nucleotides. Different prebiotic setups have successfully yielded nucleobases (4)(5)(6), sugars and their phosphate derivatives (7,8); however, the assembly of these subunits to constitute the first nucleotides has presented several obstacles (2,9).

Several strategies to produce nucleotides under prebiotic conditions have been assessed. Among them, the nucleoside phosphorylation has been one of the most widely studied. From previous works, pyrimidine and purine nucleosides were phosphorylated after a treatment with different phosphate minerals (10), or in phosphate solutions with the addition of phosphorylating agents such as urea (11), another prebiotically relevant substrate. In 2009, Powner et al. (12) proposed an elegant synthetic pathway that involved the reaction between prebiotically relevant non-canonical subunits of the nucleotides that led to pyrimidine nucleotides. Along very different synthetic routes, experiments performed by Becker et al. (13) in 2016, successfully explored purine nucleoside synthesis following the formamidopyrimidine (FaPy) pathway. A coherent prebiotic picture, encompassing the formation of both purine and pyrimidine nucleotides in similar conditions has been one of the largest current challenges in the origin of life field. In this regard, two recent works by Stairs et al. (14) and Kim and Benner (15) addressed the synthesis of both types of ribonucleotides along different chemical routes.

In current living organisms and metabolic systems, the de novo and salvage synthesis of purine and pyrimidine nucleotides implies a direct reaction between the nucleobase and a molecule of PRPP, a biological metabolite highly activated in C1' due to the presence of the pyrophosphate group. This reaction is mediated by the enzymatic catalysis of phosphoribosyltransferases specific of each nucleobase (16,[START_REF] Nyhan | Nucleotide Synthesis via Salvage Pathway[END_REF]. Primitive biochemical systems have evolved to modern ones improving their performance and selectivity: from this point of view, an interesting question is to understand whether features observed in contemporary enzymatic reactions have been conserved from prebiotic non-enzymatic pathways (chemiomimesis). In particular, one can wonder whether PRPP also played a crucial role as precursor of nucleotides in the context of a RNA world [START_REF] Kawamura | [END_REF]. In this work we examined this hypothesis (Fig. 1), which is supported by the possibility to generate PRPP as a ribose phosphorylation product (19,20), and by its plausibility to operate as a reactant toward nucleotides (20) and nucleotide mimics (ribose-pyrazole derivatives) (21). We considered mild hydrothermal conditions, that are compatible with different proposed origin-of-life scenarios, including hot surface lagoons/lakes and submarine or continental hydrothermal systems (22). We underline that we deliberately choose the simplest possible setting that allows to prove the general thermodynamic viability of the ribonucleotide synthesis reactions.

To address this crucial question, a combination of quantum-based computer simulations and experiments was carried out. Ab initio computer simulations in prebiotic chemistry under aqueous solutions (23) are very challenging due to the large system size (several hundred atoms), the bulk conditions that have to be taken into account (such as temperature, pH, ionic strength, and so forth), and especially because of the long time scale of chemical reactions. To overcome the latter obstacle, we exploited a new general approach able to simulate complex reaction mechanisms in solution, combining enhanced sampling techniques (such as metadynamics and umbrella sampling (24-26)) with topological coordinates capable of tracking changes in the chemical bond network (27). After elucidating the reaction mechanism, we could extract useful information about the environmental conditions (temperature and pH) and free energy barriers, which were further exploited to successfully test -by means of Mass Spectrometry (MS) -the formation of nucleotides in experimental conditions matching those of simulations.

Results and discussion

Ab initio simulation of nucleotide synthesis

We have assessed the formation of uridine monophosphate (UMP) and adenosine monophosphate (AMP) starting from α-PRPP and uracil or adenine, respectively, simulating easilyaccessible hydrothermal/hot lagoon conditions (T = 400 K, neutral pH), in pure water without any catalyst. Our protocol employs a preliminary exploration of the reaction pathways with metadynamics, followed by the collection of extended statistics exploiting committor analysis and umbrella sampling (see Materials and Methods).

We have discovered a simple pathway for the association of PRPP with the nucleobase, where the nitrogen forming the glycosidic bond is deprotonated. In both purine and pyrimidine cases, the reaction follows a nucleophile substitution S N 2 mechanism, with an inversion of chirality of the anomeric carbon. The free energy profile of the reaction is reported in Fig. 2, together with the detailed atomic configurations at key moments. The transition state exhibits a trigonal bipyramidal geometry, where C1' becomes planar with O and N nucleophiles from pyrophosphate and nucleobase molecules at a distance of 2.08 and 2.40 Å, respectively, for the pyrimidine, and 2.12 and 2.84 Å for the purine. The reaction connecting the deprotonated nucleobase with the nucleotide follows Hammond's postulate, i.e., the transition state is more similar to the reactants than to the products, since the forward barrier is smaller than the reverse one. When comparing the biological synthesis mediated by enzymes with the simulated prebiotic reactions, it is observed that some mechanistic details are featured in both cases: for instance the nucleobase has to be deprotonated so it can perform the nucleophilic attack on the C1' (28,29); moreover around the transition state region, the C1' adopts a planar geometry during the nucleophilic attack by the base, necessary to undergo the inversion of the chirality when forming the β-nucleotide. Overall, our results suggest that a mildly basic pH (below 9, to avoid degradation of PRPP), and an excess of the nucleobase could promote the reaction, increasing the probability of reactive collisions between a deprotonated nucleobase and the PRPP. These are the conditions we therefore adopted in our experiments.

Mass spectrometry of PRPP and nucleobase samples

Experiments were carried out in order to confirm the formation of nucleotides from nucleobases and α-PRPP, as suggested by calculations. Samples of uracil or adenine in presence of PRPP were prepared under conditions similar to the ab initio molecular dynamics simulations and subsequently analyzed by Electrospray Ionization (ESI) MS. Fig. 3 and4 show the negativemode ESI spectra of solutions containing PRPP + U, and PRPP + A, respectively, after heating at 120 • C in a closed vessel for 15 minutes. The spectra were recorded after one day aging at room temperature; reducing this waiting period did not cause significant changes in the ESI results.

PRPP is a highly activated molecule, that in water solution gives rise to degradation products well-visible in the spectra: ribose-diphosphate, 1-2 cyclic PRPP, ribose-monophosphate, cyclic ribose-monophosphate, and free phosphates. The respective amounts of these species changed with time, making it difficult to record a generally-valid PRPP reference. First, we consider the reaction of PRPP with the pyrimidine uracil. In order to determine if any UMP had been formed, we recorded the ESI spectrum of a genuine UMP solution (Fig. 3A). The spectrum is dominated by the molecular peak (UMP-1H) -at 322.9 atomic mass unit (amu). Fragmentation (MSMS, Fig. 3B) of this species yielded cyclic ribose-phosphate, and then phosphate ions, confirming the assignment (see also supplementary material, Fig. S5A-D).

Zooming in on the 323 amu region of the PRPP + U solution showed that a peak was indeed present, although with a low intensity (about 1% of the maximum intensity peak, Fig. 3D).

Since the intensity of this signal was only about three times the noise level (see Fig. S7), we carried out second fragmentation of the species at 323 amu. This yielded cyclic ribosephosphate and a small amount of phosphate ions, as in the genuine UMP reference (Fig. 3B,E), thus confirming the presence of UMP in the analyzed solution.

Clearer conclusions could be drawn for the PRPP + A solution. Here, the diagnostic peak of genuine AMP is molecular (AMP-1H) -at 345.9 amu (Fig. 4A). Indeed, a peak at this value did appear in the PRPP + A solution (Fig. 4C,D) -in this case representing about 5% of the maximum intensity peak -and its assignment was confirmed by fragmentation of the peak at 345.9 amu, Fig. 4B,E and supplementary material, Fig. S6A-C.

The relative intensity is given only as an order of magnitude, since ESI is difficult to quantify. However, it does seem that AMP is formed from PRPP + A in higher yields than UMP from PRPP + U, all other things being equal. This is consistent with the result of our simulations: the free energy difference between the products of the reaction (including the nucleotide) and the initial reagents is negative in both cases, but higher in absolute value for AMP formation, as compared to UMP formation (Fig. 2). Thus, the formation of AMP was expected to be more favored.

We applied another mass spectrometry technique, Matrix Assisted Laser Desorption/Ionization (MALDI), to check the conclusions reached by ESI. In fact, the results were not conclusive for experiments with uracil, perhaps due to sensitivity limitations. For instance, the negative-mode reference spectrum of genuine UMP has a peak at 323 amu corresponding to the singly deprotonated (UMP-1H) -anion (supplementary material, Fig. S8), but the spectrum of a heated PRPP + U solution did not show this signal (although, puzzlingly, there was a peak at 324 amu, reported on supplementary material, Fig. S9). Experiments with adenine were more conlusive. The genuine AMP reference (supplementary material, Fig. S10) in the positive mode shows two peaks at 348 and 370 amu that can be assigned to AMPH + and AMPNa + . For the PRPP + A solution, these two peaks are clearly visible in the blown-up region shown in supplementary material, Fig. S11. In addition, and quite unexpectedly, there are also significant signals that could be attributed to species bearing more phosphoryl groups: ADPH + (or one of its isomers) at 428 amu, with ADPNa + at 450 amu and (ADPNa 2 -1H) + at 472 amu, and ATPH + at 508 amu, with ATPNa + at 530 amu (see supplementary material, Fig. S12-S13). These data were confirmed by negative-mode MALDI where (AMP-1H) -was detected at 346 amu, as well as (ADP-1H) -at 426 amu and (ATP-1H) -at 506 amu as reported in supporting information (see supplementary material, Fig. S14-S15).

Following glycosylation, nucleobases could in principle bind to ribose through other sites than N9 for purines and N1 for pyrimidines. We note, however, that in deprotonated uracil the most stable resonance structure features the negative charge on N1 (30), indicating that the glycosidic bond is more likely to be formed on that site. In the case of adenine, a specific regioselectivity study should be performed, to distinguish the different possible sites for the formation of the glycosidic bond.

In summary, mass spectroscopic data show the formation of purine and pyrimidine nucleotides in PRPP + nucleobase hydrothermal solutions. In the case of PRPP + U, although MALDI, likely because of sensitivity reasons, does not detect the UMP nucleotide, ESI strongly indicates its formation. We explicitly verified using 31 P NMR that the ESI-detected UMP (as well as AMP) was not present as a contaminant in the PRPP reagent employed for the synthesis reactions (see supplementary material, Fig. S16). In the case of PRPP + A, both ESI and MALDI clearly detect the AMP nucleotide, probably with larger yields than UMP, consistently with free-energy calculations. Furthermore, ab initio simulations predict that the largest barrier for nucleotide synthesis (Fig. 2) is of about 22 k B T (T = 400 K), pointing to a relatively fast kinetics that is compatible with the time scale of nucleotide synthesis in our experiments. The relatively low yields seem to be due mainly to the competition between nucleotide formation and PRPP degradation reactions, by hydrolysis and cyclic phosphate formation. Further work is required to better understand this competition and determine what settings favor most nucleotide formation. For instance, in plausible prebiotic conditions, mineral surfaces are known to play a crucial role in favoring both the synthesis and stability of PRPP (20), while wet-dry (hydrationdehydration) cycles could favor nucleotide accumulation and further polymerization (31).

Conclusions

In conclusion, we present hereby two main discoveries. First, we identify a new synthetic route spontaneously leading to both purine and pyrimidine ribonucleotides, starting from precursors of both prebiotic and biological relevance. This route is compatible with several hydrothermal environments, from submarine vents to surface lagoons and lakes, most likely present in the early Earth as well as on other planets. In this context, mineral surfaces could have played an important role in the synthesis and stabilization of PRPP, as well as in improving nucleotide condensation yields by reducing water activity (20).

Second, we demonstrate, with atomic-level precision, that this spontaneous prebiotic synthetic pathway is very similar to the biological one, to the point of sharing the same precursors, the same steps, and the same chirality-inverting character (Fig. 1). This correspondence is highly suggestive of an evolution process towards complex life forms conserving the basic features of prebiotic ribonucleotide synthesis. Building upon these findings, further studies will be able to assess under which conditions the reactions can proceed up the ladder of complexity, starting from the oligomerization of nucleotides, until the possible inception of a fully-functional RNA world.

Materials and Methods

Born-Oppenheimer Molecular Dynamics (BOMD)

Ab initio calculations were based on the Density Functional Theory with the exchange-correlation Perdew-Burke-Ernzerhof functional (32), including Grimme's dispersion corrections (33) and the Martins-Troullier pseudopotentials (34) for C, N, O, Na, P and H atoms, as implemented in the code CPMD 4.1 [START_REF] Corp | Max Planck Institute Stuttgart[END_REF]. BOMD simulations were performed in the canonical ensemble (NVT) at a temperature of 400 K, with a time step of ∼0.48 fs, in a periodically repeated cubic box of approximately 4.296 and 4.292 nm 3 , and a density of 0.9888 and 0.9984 g/cm 3 for UMP and AMP respectively. Temperature was controlled by the Nosé-Hoover thermostat with a chain length of four and with a frequency of 3000 cm -1 . The electronic convergence criterion has been set to 10 -5 a.u. based on benchmark NVE simulations. The simulation boxes comprise a PRPP molecule with charge 5-and a uracil/adenine molecule (reactants) or a nucleoside monophosphate molecule (UMP or AMP) and a pyrophosphate ion in its HP 2 O 3- 7 form (products), with additional sodium counterions and 108 explicit water molecules, for a total of 371 atoms in the pyrimidine setup and 374 for the purine.

Enhanced Sampling

BOMD simulations were coupled with enhanced sampling techniques. In the first place we performed metadynamics simulations (24), employing as set of reaction coordinates the path collective variables (pathCV) s and z by Branduardi et al [START_REF] Branduardi | [END_REF]. In general, s indicates the progress along a putative pathway composed by a discrete sequence of atomic configurations, while z represents the distance from the putative pathway (complete definitions are provided in supplementary material). In this work, we employed only two reference configurations, corresponding to reactants and products, in this way avoiding to bias the observed reaction mechanisms with hypotheses on the path. In addition, we used a topological metric particularly suited to chemical reactions in solution (27), accurately tracking changes in the chemical bond network passing from reactants, through intermediates, until the products. This metric comparing structure R(t)

with reference structure R k is defined as D(R, R k ) = IS C IS -C k IS 2
, where coordination numbers of a specific atom I with the atoms of species S appear.

In the second stage, we carried out a committor analysisto identify the transition state of each reaction and to get more insights about its mechanism. We extracted a set of about 10 atomic configurations from reactive metadynamics trajectories, and we performed up to 20 independent unbiased BOMD trajectories of 0.5 ps starting from each configuration, verifying whether they fell into reactants or products basins. A transition state configuration is identified as committed to both basins with a probability of 50±10%.

In the final stage, a series of umbrella sampling (25) simulations were performed, systematically restraining the s coordinate at different locations along the reaction pathway obtained by committor analysis using an harmonic potential (all details are reported in supplementary material, Table S5). The trajectory length of each window was 10 ps, for a total simulation time of 250 and 220 ps for the pyrimidine and purine cases, respectively. The weighted histogram analysis method (37) finally yielded the free energy profile of the reactions. All enhanced sampling 3. Mass spectrometry (ESI-IT-MS -negative mode) of (A) UMP reference, with the peak at 322.9 amu corresponding to (UMP-1H) -, (B) fragmentation of peak 322.9 amu in UMP reference, (C) complete spectrum for the sample PRPP + uracil, (D) zoomed spectra around UMP region (317-330 amu), and (E) fragmentation of peak at 322.9 amu from sample PRPP + uracil. The detailed labeling of peaks is explained in Supporting Information. Fig. 4. Mass spectrometry (ESI-IT-MS -negative mode) of (A) AMP reference, with the peak at 345.9 amu corresponding to (AMP-1H) -, (B) fragmentation of peak 345.9 amu in AMP reference, (C) complete spectrum for the sample PRPP + adenine, (D) zoomed spectra around AMP region (337-358 amu), and (E) fragmentation of peak at 345.9 amu from sample PRPP + adenine. The detailed labeling of peaks is explained in Supporting Information.

FigFig. 1 .

 1 Fig S11. MALDI positive mode PRPP + A, formation of AAP Fig S12. MALDI positive mode PRPP + A, formation of ADP Fig S13. MALDI positive mode PRPP + A, formation of ATP Fig S14. MALDI negative mode PRPP + A, zoomed region 304-397 amu Fig S15. MALDI negative mode PRPP + A, zoomed region 400-520 amu Fig S16. 31 P NMR for genuine AMP, UMP and PRPP Other supplementary material: Movie S1. reaction PRPP + U Other supplementary material: Movie S2. reaction PRPP + A

Fig. 2 .

 2 Fig. 2. Relative free-energy values for reactants (A), intermediates (B), transition states (C), and products for the nucleotide synthesis (D). Free-energy levels for uracil (pyrimidine) and adenine (purine) pathways are depicted in green and orange respectively. "N*" represents the nitrogen from the nucleobase participating in the glycosidic bond. Relative free-energy values for reactants (A), are estimated based on the pKa values of the nucleobases (39, 40). Structural representation of reactants (B U and B A ), transition states (C U and C A ), and products (D U and D A ) during the nucleotide formation from PRPP and uracil or adenine, respectively. The transition states display the typical geometry of S N 2 mechanisms.

Fig.

  Fig.3. Mass spectrometry (ESI-IT-MS -negative mode) of (A) UMP reference, with the peak at 322.9 amu corresponding to (UMP-1H) -, (B) fragmentation of peak 322.9 amu in UMP reference, (C) complete spectrum for the sample PRPP + uracil, (D) zoomed spectra around UMP region (317-330 amu), and (E) fragmentation of peak at 322.9 amu from sample PRPP + uracil. The detailed labeling of peaks is explained in Supporting Information.

Acknowledgements Funding: This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02, within the framework of the cluster of excellence MATériaux Interfaces Surfaces Environnement (MATISSE) led by Sorbonne Université. We acknowledge calculations performed on the Ada cluster at IDRIS, (Orsay, France) under GENCI allocations 2015-091387 and 2016-091387; and Curie Supercomputer at TGCC under GENCI allocations t201609s042 and A0010910143.

Sample preparation

All data presented in this work were recorded with controlled and pure PRPP purchased from Sigma Aldrich. We further verified the lack of any nucleotide contaminant by 31 P NMR (see supplementary material, Fig. S16). 0.1 mol/L solutions of PRPP, UMP, and AMP were prepared in pure water as reference. Mixtures of PRPP and bases (uracil or adenine) were prepared by dissolving 12.5 mg of PRPP (0.05 mol/L) with 11.5 mg of uracil (0.2 mol/L) or 13.5 mg of adenine (0.2 mol/L) in 500 µL of pure water. When needed, the pH was adjusted to 8.5 by controlled addition of NaOH. Before mass spectrometry analysis, solutions are filtered with 2µm paper. The reactions are carried on in a glass vial with milliQ water. Then, the reactants are thermally activated at 120 • C during 15 minutes.

ElectroSpray Ionization -Ion Trap -Mass Spectrometry (ESI-IT-MS)

Samples were diluted to the micromolar range in a solution of 50% acetonitrile, 49.8% water and analyzed by direct infusion in a HCT Ultra PTM Discovery mass spectrometer from Bruker Daltonics (Germany) equipped with an electrospray ion source. The nebulizer gas pressure was 11 psi and the spray voltage was ∼3.5 kV. The drying gas flow was 5 L/min and the temperature was 300 • C. Spectra were acquired in negative ion MS mode over a 50-2000 m/z range until the ion charge control target had reached 90000. The MS acquisitions were carried out at enhanced resolution (8100 m/z units per second). External calibration was performed with the ESI-T Tuning Mix (Agilent Technologies). The instrument was controlled using EsquireControl 6.2 software (Bruker Daltonics) and mass spectra were processed using DataAnalysis 4.0 SP2 software (Bruker Daltonics).

MALDI-TOF/TOF

Mass spectra were generated using a 4700 Proteomic Analyzer MALDI-TOF/TOF (Applied Biosystems) instrument fitted with a Nd:YAG laser (λ = 355 nm; pulse duration = 4 ns; repetition rate = 200 Hz). A dilute solution was vortexed for 1 min at high speed prior to deposition on the MALDI plates. All MALDI-TOF spectra, resulting from the average of a few tens to a hundred laser shoots, were obtained in positive and negative ions using the reflector mode in the m/z range of .

Supplementary Materials

Supplementary Text: 1.1 Enhanced sampling Table S1. Coordination matrix for PRPP + U solution Table S2. Coordination matrix for UMP + PPi Table S3. Coordination matrix for PRPP + A solution Table S4. Coordination matrix for AMP + PPi Table S5. Umbrella sampling parameters Supplementary Text: 1.2 Assignment of MS peaks Compelling interests: The authors declare no competing financial interests.

Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.


Laurent Nony
email: laurent.nony@im2np.fr.

Franck Bocquet

Christian Loppacher

Thilo Glatzel

On the relevance of the atomic-scale contact potential difference by Amplitude modulation-and Frequency modulation-Kelvin probe force microscopy

The influence of short-range electrostatic forces on the measured local Contact Potential Difference (CPD) by means of Amplitude Modulation-and Frequency Modulation-Kelvin Probe Force Microscopy (AM-and FM-KPFM) is discussed on the base of numeric and analytic descriptions of both methods. The goal of this work is to help interpreting recent experimental results reporting atomically-resolved CPD images, in particular on bulk insulating samples. The discussion is carried out on the base of spectroscopic curves. The expression of the bias-dependent electrostatic force derives from a previous work and is estimated between a tip with simple geometry and the (001) facet of a perfect alkali halide single crystal. The force, with a short-range character, scales as a second-order polynomial function of the bias voltage. It is stated that the linear term is responsible for the occurrence of the atomic-scale CPD contrast, while the quadratic one, involving the sample polarisation, accounts for the detected signal by the KPFM methods. Nevertheless, analytic and numeric approaches stress the influence of the linear term on the measured CPD which intrinsically hinders the possibility to perform quantitative CPD measurements, but also makes the measured "pseudo-CPD" strongly deviating from the surface potential. Hence, in the short-range regime, AM-or FM-KPFM measurements neither reflect the CPD nor the local surface potential, but rather an effective value which is convoluted by the geometric parameters of the tip, the so-called local CPD. It is also stated that the local CPD measured by means of AM-or FM-KPFM differs when sub-nanometer vibration amplitudes of the cantilever are used. Otherwise, AM-and FM-KPFM measurements should be almost similar. At last, the influence of long-range, capacitive, electrostatic forces is discussed in conjunction with the short-range ones. This allows us to draw conclusions regarding the distance dependence of the local CPD which then exhibits a resonant behavior as a function of the tip-surface separation. This phenomenon is expected to play a role in the KPFM imaging process.

I. INTRODUCTION

Kelvin Probe Force Microscopy (KPFM) is a scanning probe method 1 which is based on the detection and the dynamic compensation of the electrostatic forces occurring between a nano-tip and a sample when they are electrically connected. These forces stem from intrinsic work function differences between the tip and the sample surface 2,3 , capacitive contributions, as well as from trapped charges at the tip apex and/or the surface upon preparation procedures 4 . The KPFM method actually provides access to the contact potential difference (CPD) 5 , i.e. the work function of the sample referenced to the work function of the tip: qV cpd = ∆φ = φ tip -φ sample , q being the elementary electrical charge.

The technique is based on the regular noncontact-Atomic Force Microscopy (nc-AFM) operating mode, which is strongly sensitive to electrostatic forces. Two KPFM-based techniques provide facilities to map the spatial variations of the CPD on the nanometer scale, namely Amplitude Modulation-KPFM (AM-KPFM, refs. [6,7]) or Frequency Modulation-KPFM (FM-KPFM, ref. [8]). The connection between KPFM and nc-AFM has brought a complimentary information to the usual topography and more controversial dissipation (socalled damping) channels, and has already been used to map the chemical identity of surface atoms 9 . KPFM early proved its ability to map the spatial variations of the CPD on the nanometer scale with a resolution of few mV 6,10,11 . Different groups reported atomicallyresolved CPD images, some including even the influence of atomic adsorbates 8,9,12,13,14,15 .

More recently, the CPD atomic-scale contrast was reported on metallic and ionic thin films on InSb 16 , on oxides 17,18 and on ionic bulk insulating surfaces 19 . The former results suggest that, very close to the surface, the CPD rather has a local character 20 . Among the former references however, despite the consistency between the lateral periodicity of the local CPD images and the surface lattice constant, the values rarely fit to the CPD values reported by other experimental methods or with the theoretical predictions 8,12,14 . For instance, on the complex Si(111)5 √ 3 × 5 √ 3-Sb surface, Okamoto et al. concluded that the CPD images mainly reflect the electrostatic force distribution on the surface rather than the work functions distribution 9,15 . It has also been mentioned that the measured CPD was dependent on the tip-surface distance 15,16,21 . Thus, a cross-talk between the topography and CPD images is very likely, as soon as the topography feedback circuit is engaged.

These results make the interpretation of the CPD atomic-scale contrast difficult to explain and the measurements ultimately questionable. Besides, the CPD being intrinsically a macroscopic property of the tip-surface system, there is currently an intense debate dealing with the apparent contradiction between CPD and local CPD measurements, which actually raises the questions of the origin and the reliability of atomic-scale resolution in KPFM 16,18,19 .

The difficulty for interpreting the experimental results stems from the lack of theoretical description for the short-range electrostatic forces which are intricately involved in the atomic-scale contrast and the way these forces are processed by the KPFM control electronics. In order to address these questions, we first have focussed on the CPD atomic-scale contrast reported experimentally on the KBr(001) insulating surface 19 . In this preliminary work, an analytic model for the short-range electrostatic interaction between a biased tip and the KBr surface was developed.

In the present work, the influence of such a short-range electrostatic force on the measured local CPD is discussed by means of a numerical implementation of the AM-KPFM and FM-KPFM setups within our nc-AFM simulator 22 . The numerical results are completed with an analytical approach providing approximated expressions of the modulated components used by the AM-and FM-KPFM methods as a function of the DC bias voltage. The following section reminds the AM-and FM-KPFM setups, their operational mode and their experimental implementation, which were accurately duplicated on the numerical level. In section III, the analytical model for the short-range electrostatic force is briefly introduced and we stress the influence of the dynamic polarisation of the ionic crystal, which is for a large part responsible for the occurrence of the atomic contrast. In sections IV A and IV B, the analytic and numeric approaches are discussed in respect with the CPD atomicscale contrast. In section V, we extend the framework of the atomic-scale CPD contrast to long-range electrostatic forces. This allows us to draw conclusions regarding the distance dependence of the local CPD, which is not possible when considering uniquely short-range electrostatic forces.

II. AM-AND FM-KPFM EXPERIMENTAL METHODS

KPFM relies on the dynamic compensation of the electrostatic forces arising between the tip and the sample, hereafter referred to as the electrodes. For that purpose, an amplitude modulation technique is used. The electrostatic force is triggered when modulating the bias voltage which is applied between the electrodes:

V b (t) = V DC + V mod cos(2πf mod t) (1)
When dealing with metallic electrodes and for large tip-surface distances compared to atomic dimensions, the electrostatic part of the total force can be written as:

F es (t) = - 1 2
∂C(z) ∂z (V b (t) ± V cpd ) 2 , (2)
C(z) being the expression of the tip-surface capacitor, influenced by the geometry of both electrodes 23,24 . The ± sign depends upon the bias voltage is applied to the tip or to the sample. When applied to the tip (sample grounded), the measured DC potential is the opposite of the CPD, hence one should write V b (t) + V cpd , otherwise the opposite (tip grounded). In the following, we will assume that the tip is grounded and will write V b (t) -V cpd .

The above expression of the force gives three temporal components, the first one being static, the second and third ones being amplitude-modulated at f mod and 2f mod , respectively.

The relevant component for the KPFM technique is the one modulated at f mod :

F f mod es (t) = ∂C(z) ∂z (V DC -V cpd )V mod cos(2πf mod t) (3)
This component is fed into a dual phase lock-in amplifier (LIA) which detects the amplitude of the signal, A K , from the measurement of the in-phase and out-of-phase components X and Y , respectively. The signal X is then injected into a proportional/integral controller, the output of which supplies the proper DC bias voltage to minimize, or ideally cancel X and hence A K , i.e. also F f mod es (cf. equ.3). It is important to notice that X is supposed to be used instead of A K in the Kelvin controller, because it can become negative and thereby handle negative error signals, while A K not. Thus, according to equ.3, the output of the controller is the DC bias which compensates the CPD:

V (c)
DC = +V cpd . The experimental implementation can be done in two ways, both of them relying on the detection and subsequent cancellation of A K , but differing regarding the signal that carries the modulation. In FM-KPFM, the modulation is performed at a few kHz and detected in the frequency shift signal (∆f ) provided by the phase-locked-loop 25 (PLL). In AM-KPFM, the modulation frequency matches the second flexural eigenmode of the cantilever 7,26 , with typical frequency 27 : f 1 ≃ 6.24f 0 . The experimental setups have been depicted in fig. 1. They are based on the well-established nc-AFM operating mode which relies on the PLL-controlled excitation of the cantilever 22,28 .

Equs.2 and 3 show that the KPFM methods intricately relies on electrostatic forces scaling quadratically with the bias voltage, otherwise no CPD contrast would be detected by the experimental method. Besides, atomically-resolved CPD images testify that these forces have a short-range character. However, equ.2 is improper to describe this contrast, first because the distance dependence, i.e. the expression of the tip-surface capacitor, has usually not a short-range character and second because this expression derives from a classical electrostatic approach between continuous-like metallic bodies and hence does not take into account atomic fluctuations of the surface potential of the insulating sample. Another fundamental aspect of the problem stems from the influence of the technique onto the measured CPD. In our previous work 19 as well as in the work by Krok et al. 16 , there are strong evidences that short-range electrostatic forces scale not only quadratically, but also linearly with the applied bias voltage. Assuming a linear term in equ.2:

F es = C 2 (V b (t) -V cpd ) 2 + C 1 (V b (t) -V cpd ), then the modulated component of the force is modified consistently: F f mod es = [C 1 + 2C 2 (V DC -V cpd )]V mod cos(2πf mod t). The compensated CPD becomes: V (c) DC = V cpd -C 1 /(2C 2 )
, which fakes its physical content. The physics that it carries now differs from the actual CPD and depends on the coefficients C 1 and C 2 . Obviously, this remains true when other power laws of the bias are to be considered in the expression of the electrostatic force.

III. SHORT-RANGE ELECTROSTATIC FORCE MODEL

In ref. [19], Bocquet et al. proposed an analytical approach to the short-range electrostatic force between a nanoscopic tip and the (001) surface of a perfect alkali halide single crystal of KBr, i.e. not including residual net charges. The sample is several millimeters high, which places the metallic counter-electrode far from the tip in comparison with other setups. The two former assumptions insure that no Coulombic nor capacitive long-range electrostatic forces act onto the tip in the process that is being discussed. The tip consists of two embedded half-spheres with different radii. The bigger radius R stands for the mesoscopic body of the tip (R ≃ 5 nm). The sphere with radius R a ≪ R, which is half-embedded within the macroscopic body, rather stands for a microscoscopic nanoasperity supposed to FIG. 1: Scheme of the combined nc-AFM / KPFM experimental setup which was duplicated numerically. AM-or FM-KPFM mode is selected upon the position of the switch. The PLLcontrolled excitation scheme was detailed in ref. [22].

favor the occurrence of the atomic-scale contrast (R a ≃ 1 Å, cf. fig. 2a). With this simple geometry, the expression of the force can be split into two major contributions:

F es (V b , z) = F (1) mµ + F (2) µ = C 1 (V b (t) -V cpd )e -αz + [C 0 + C 2 (V b (t) -V cpd ) 2 ]e -2αz , (4)
α is a coefficient which scales with the lattice constant of the crystal, a, according to α = 2 √ 2π/a. The coefficients C 0 , C 1 and C 2 are written:

C 0 = - R 2 q 2 2ǫ 0 a ′4 A (2) e -2αRa and C 2 = 2C 0 2πχ d a ′ qR 2 (5)
and:

C 1 = - 3R 2 a q ǫ d a ′2 R D (1) e -αRa [cos( x 0 ) + cos( y 0 )] (6)
R and R a are the tip's radii defined above. ǫ 0 , ǫ d and χ d are the vacuum dielectric permittivity, KBr effective dielectric permittivity and KBr dielectric susceptibility, respectively. In the former reference, we had set: ǫ d = 4.38 and χ d = 9 10 -39 Fm 2 . A (2) and D (1) are two numeric coefficients depending on the values of radii R and R a in a complex manner (cf.

equs.A2 and A6 in ref. [19]). x 0 = 2πx 0 /a ′ and y 0 = 2πy 0 /a ′ are the reduced coordinates of the position of the tip onto the KBr(001) surface. Setting x 0 = y 0 = 0 or x 0 = y 0 = a ′ /2 locates the tip on top of an anion or of a cation, respectively. Thus, the coefficient C 1 carries the lateral modulation of the force over the crystal surface.

F (1)
mµ depicts the coupling between the microscopic structure of the tip apex and the capacitor consisting of the tip, the ionic crystal and the counter-electrode. It originates from the influence of the KBr Madelung surface potential (MSP) onto the surface charge density that develops within the microscopic nanoasperity due to the bias voltage V b . F

(2) µ , depicts the influence of the MSP onto the mesoscopic part of the tip, independently from its microscopic structure. The exponential distance dependence, which stems from the MSP, shows that the force has a short-range character.

F (1)
mµ has the lateral periodicity of the MSP (cf. coefficient C 1 ), whereas

F (2)
µ exhibits no spatial dependence and merely acts as a static component, which shifts the total force. Therefore the CPD atomic-scale contrast is mainly carried by the term F

(1) mµ . The total force is in the range of ten pico-Newtons. Let us finally note that the quadratic bias voltage dependence of F

(2) µ occurs if and only if the ionic polarisation of the sample due to the influence of the tip/counter-electrode capacitor is explicitly considered. By ionic polarisation, it is meant a net vertical displacement of the ions from their equilibrium positions within the crystal upon their sign, which changes the strength of the MSP 29 . This displacement follows the bias modulation dynamically.

The above elements allow us to draw two major conclusions: (1)-the measurement of the local CPD by means of KPFM is partly made possible via the ionic polarisation of the sample (required quadratic bias voltage dependence) and ( 2)-the occurrence of the atomicscale contrast, indeed consistent with the ionic positions, relies on a linear bias voltage dependence of the short-range electrostatic force. Therefore, according to the discussion of section II, it is expected that the compensated CPD, either measured by AM-or FM-KPFM, does not match the actual tip-surface CPD.

FIG. 2: a-Schemes of the three tips used for the analytic and numeric calculations, properly scaled. The set of parameters (R, R a ) in Å are: tip1(50, 1); tip2(120, 1.5); tip3(50, 1.1). For the analytic calculations, the sets of coefficients (A (2) , D (1) ) are: tip1(-2.5,-15); tip2(-1.1,-25); tip3(-2.5,-20). b-Scheme of the vibration of the cantilever at the 1 st and 2 nd eigenmodes, corresponding to instantaneous positions z 0 (t) and z 1 (t), respectively. The instantaneous tip-surface separation

is z(t) = D -z 0 (t) -z 1 (t).

IV. ANALYTIC AND NUMERIC APPROACHES

A. AM-KPFM

Analytic approach

In ref. [19], a self-consistent approximated expression of the vibration amplitude of the second flexural eigenmode of the cantilever was proposed. Let z 1 (t) be the instantaneous deflection of this mode (cf. fig. 2b) and A 1 its vibration amplitude. It is ruled by the nonlinear equation of motion:

z1 (t) + ω 1 Q 1 ż1 (t) + ω 2 1 z 1 (t) = F ext m 1 + F es (V b , z) m 1 , (7)
where Q 1 , ω 1 = 2πf 1 and m 1 are the quality factor, the resonance pulsation and the effective mass of the mode, respectively. F ext stands for the actuation force of the fundamental flexural mode of the cantilever (resonance frequency f 0 ). The resonance amplitude of this mode is A 0 . Owing to the large difference between f 0 and f 1 , F ext does not influence z 1 (t). Let z 0 (t) be the instantaneous position of the fundamental mode. Thus, the instantaneous tipsurface separation to be considered for F es (V b , z) is: z(t) = D -z 0 (t) -z 1 (t), D being the separation between the surface and the equilibrium position of the cantilever at rest (cf.

fig. 2b). Owing to the exponential distance-dependence of F es (V b , z), the second mode is not actuated harmonically, but rather stepwise-like upon z(t). This makes the structure of the actuation of the second mode rather complex, but nevertheless enough to permit the development of the steady state, which merely stems from the large value of Q 1 (a few thousands).

When linearizing equ.7 with respect to z 1 (t) and assuming it as harmonic, an approximate expression for A 1 can be derived (cf. ref. [19] for details):

A 1 = l 0 (l 3 -l 2 ) 2 + l 2 1 l 2 1 + l 2 2 -l 2 3 ( 8
)
Except l 0 , coefficients l 1 , l 2 and l 3 are reported in the appendix. l 0 is written:

l 0 = C 1 a 0 m 1 e -αz min V mod + 2 C 2 b 0 m 1 e -2αz min (V DC -V cpd )V mod (9)
z min = D -A 0 stands for the minimum tip-surface separation, i.e. at the lower turning point of the tip when neglecting A 1 with respect to A 0 , which is practically always correct (see below). a 0 and b 0 are the the zero-order Fourier coefficients of the functions exp{-αA 0 [1cos(2πf 0 t)]} and exp{-2αA 0 [1cos(2πf 0 t)]} that occur in the expression of the force when expanding z(t). The n th Fourier coefficients are: a n = 2 exp(-αA 0 )I(n, αA 0 ) and b n = 2 exp(-2αA 0 )I(n, 2αA 0 ), I being the modified function of Bessel of the first kind 19 . Hence, the condition on the bias to nullify A 1 , i.e. the modulated component at f mod in AM-KPFM, is simply l 0 = 0, which yields:

V (c) DC = V cpd - C 1 2C 2 a 0 b 0 e +αz min (10)
The above expression shows that the compensated CPD has the lateral periodicity of the MSP (coefficient C 1 ). However, it also depends on the tip parameters R and R a . This result shows that the experimental measurements of the local CPD by AM-KPFM neither reflects the actual tip-surface CPD, nor the value of the local surface potential, but rather an effective value that is convoluted by the tip geometry. This originates from the linear bias dependence in the expression of the force, as stated before. Nevertheless, the lateral periodicity of the MSP is preserved in the CPD contrast, which makes CPD differential measurements among atomic sites reliable.

Equ.10 shows that, above an anion, the compensated CPD becomes more and more negative as the tip-surface separation is increased 30 . This rather surprising behavior stems from the faster distance dependence of F

(2)

µ compared to F (1)
mµ , as discussed in ref. [19]. Therefore, the distance dependence of the compensated local CPD depends much on the tip geometry.

Obviously, this behavior cannot stand for any tip-surface separation owing to the exponential decay of the short-range electrostatic force which triggers the motion of the cantilever's 2nd eigenmode and hereby, the detection of the local CPD. If the force is too weak, no CPD can be detected anymore. However, within a narrow range of tip-surface separations, typically a few angströms, it might be detectable. Nevertheless, the above discussion is not complete because it is insufficient to consider the short-range electrostatic force as the main triggering source of the cantilever since it is known that long-range capacitive forces always occur on the experimental level. In section V, the distance dependence of the local CPD is discussed in more detail within the framework of short-and long-range electrostatic interactions.

Numeric approach

The numerical implementation of the AM-KPFM setup has been performed within the code of the nc-AFM simulator. In particular, although the resonance frequency of the second flexural mode is in the MHz range (f 1 = 6.24f 0 = 6.24 × 150 kHz ≃ 940 kHz), we have kept the sampling frequency standing for the analog parts of the electronics constant, namely f s 2 = 400 MHz (cf. ref. [22] for details). This is still sufficient to integrate the differential equation of motion of the mode with an error kept low enough. The implementation sticks to the experimental setup shown in fig. 1. We have used a simple first-order, high-pass filter with a 200 kHz cutoff frequency. The dual-phase LIA has a 10 kHz bandwidth. The numerical implementation is similar to the one used in the nc-AFM simulator for monitoring the phase lag between the excitation and the cantilever's response. It provides

A K = 2 √ X 2 + Y 2 ,
which in the case of the AM-KPFM matches A 1 , the vibration amplitude of the second eigenmode. The Kelvin controller has a 2.5 kHz bandwidth. It is a standard proportional and integral controller, the implementation of which is similar to the distance and amplitude controllers. In the present work, we merely have focussed at A 1 vs. V DC curves, so-called spectroscopic curves. Therefore, the measurements are performed at a single (x, y) position and the Kelvin controller is not engaged. The determination of the influence of the Kelvin controller to the measured local CPD will be addressed in a future work.

At each time step, the code integrates the differential equation for the second flexural mode (equ.7) in parallel to the equation of motion for the fundamental flexural eigenmode, which has the same form. The PLL-excitation scheme ensures that the latter mode is continuously actuated at its resonance frequency which shifts as the tip is brought closer to the surface. According to equ.4, the two equations of motion are nonlinear and coupled via the instantaneous separation between the tip and the surface z(t) defined before.

In order to focus on the origin of the local CPD, no additional forces like long-range Van der Waals or electrostatic ones, or short-range chemical ones, have been included in the simulations. Hence, the cantilever dynamics is merely influenced by the electrostatic shortrange force field. The force field is implemented as a 2D lookup table with the z variations sampled every 2 pm and bias voltage variations every 10 mV, the boundary values being adjustable upon needs. The sequence of simulation of an AM-KPFM spectroscopic curve is the following: the cantilever steady state is calculated at a tip-surface separation for which the interaction is zero and with V DC = 0 V. The PLL and amplitude controller are then engaged and their parameters set optimally 22 . Then, the approach is performed down to an arbitrary value of z min . During the approach, the second mode undergoes a frequency shift owing to the influence of the nonlinear interaction, like the first one. But it is important to ensure the on-resonance excitation of this eigenmode too. This requires to tune precisely the modulation frequency of the bias to recover the on-resonance excitation. Finally, the A 1 vs. V DC spectroscopic curve is acquired by continuously sweeping the DC part of the bias, first from 0 down to negative values and then upwards.

Results and discussion

Analytic and numeric curves are reported in fig. anion of the KBr surface. At similar height, the MSP of the (001) KBr facet is -440 mV. The choice of that value of z min stems from the exponential decay of the electrostatic short-range force. One has to be close enough to the surface to get a force able to trigger the mode. For that distance, the force yields -52 pN, -35 pN and -46 pN with tips 1, 2 and 3, respectively. This is also why we have set a relatively large value for the bias modulation, although this value can be significantly lowered experimentally owing to longer-range electrostatic forces. For the sake of clarity we have set the CPD of the tip-surface system equal to zero, V cpd = 0 V. The DC part of the bias is swept from -5 V to +5 V within 30 s. The values of the local CPDs are reported in table I.

Analytic and numeric curves give qualitatively similar trends, both pointing out the strong influence of the tip geometry on the measured local CPD. Although the asperity radius between tips 1 and 3 merely differ by 10% (R a = 1 and 1.1 Å, respectively), their relative local CPDs are shifted by several volts (1.5 V with the analytic approach and about 3 V for the numeric approach). This stems from the term F

(2) µ which, in the case of tip 3, is decreased because the mesoscopic part of the tip is 10 pm farther away from the surface than for tip 1. Therefore the influence of the MSP, which decays exponentially fast, onto the tip is weaker. The strength of F

(2) µ being governed by the coefficient C 2 , the above argumentation can be interpreted as if C 2 was strongly lowered between tip 1 and tip 3

(given a value of z min ). According to equ.10, this must result in a consistent change of the local CPD. Regarding tip 2, the effect is stronger even as the asperity is bigger than for the two former tips. Furthermore, the mesoscopic part of the tip being bigger, the tip-surface capacitor is larger, which strengthens the influence of the term F

(1) mµ , i.e. coefficient C 1 in equ.10. The influence of tip's radii R and R a on the measured local CPD is discussed in more detail in section V.

A straightforward consequence is that the measured local CPDs neither match the CPD of the tip-surface system (0 volt), nor they are comparable to the MSP at similar distance.

The values are systematically much over-estimated (absolute values). This was predicted by the analytical approach and is indeed confirmed numerically. Thus, it is an effect to be considered for interpreting the experiments.

The above elements must not hide the quantitative discrepancies between analytic and numeric curves. First, the analytic approach lacks in providing a value of the local CPD similar to the numerical one and second, the curvature of the curves at large bias differs between the two approaches. Regarding the numeric approach, the differential equations are integrated consistently without approximations. Therefore, we believe that the numeric local CPDs are somewhat reliable. On the contrary, the analytic approach relies on numerous assumptions, some of which being strong (harmonic behavior of the mode). This is why it is possible to derive a simple expression for V (c) DC wherein all the assumptions are actually gathered within coefficients a 0 and b 0 , i.e. the zero-order Fourier components of the short-range electrostatic force. When tuning them by only few percents, the value of

V (c)
DC is significantly modified. Furthermore, due to the linearization procedure, the analytic solution for A 1 is basically valid for a small force, i.e. in a small region around the CPD value. Therefore, the curvature of the numerical curves is likely connected to the substantial increase of the electrostatic force which modifies the excitation, and consequently the dynamics, of the second eigenmode.

B. FM-KPFM

Analytic approach

In FM-KPFM, the bias voltage modulates the electrostatic force and thus, the oscillation of the cantilever is frequency-modulated at f ′ 0 + f mod and f ′ 0 + 2f mod , as shown in the work by Zerweck et al. 25 . f ′ 0 is the resonance frequency of the cantilever close to the surface, i.e. shifted from its natural value f 0 due to the influence of other forces like Van der Waals and chemical short-range ones, but also the DC part of the electrostatic force. This means that, for a given value of z min , the frequency shift supplied by the PLL is not constant anymore, but modulated at f mod and 2f mod . It is reminded that in FM-KPFM, f mod is of about 1 kHz and hereby far from any resonance of the cantilever. Therefore, it is also very low compared to f ′ 0 and f 0 . In order to derive an analytic expression for ∆f , we use the formula introduced by Giessibl 31 . When applied to FM-KPFM, this approximated expression of ∆f is very good because, owing to the range of the electrostatic force, the perturbation of the harmonic motion of the first eigenmode of the cantilever is weak (particularly the one due to the AC part), but also because the control electronics is designed to maintain the harmonic motion of the cantilever. We then perform the calculation of ∆f by using this formalism coupled to the following quasi-static approximation: because f mod ≪ f 0 , it is assumed that the AC part of the electrostatic force is constant during one oscillation cycle of the cantilever. In ref. [16], Krok et al. have used a similar approach with a different expression for the electrostatic force, however scaling similarly with V b . Hence:

∆f = f 0 2k 0 A 0 e -αz min { C 1 (V DC -V cpd )a 1 + b 1 e -αz min [C 0 + C 2 ( V 2 mod 2 + (V DC -V cpd ) 2 )] + [C 1 a 1 + 2(V DC -V cpd )C 2 b 1 e -αz min )] V mod cos(2πf mod t)+ [C 2 b 1 e -αz min ] V 2 mod 2 cos(4πf mod t) (11)
a 1 and b 1 are now the f 0 -Fourier components (first components) of the exponential functions discussed above. The first term in equ.11 is static and gives, for V mod = 0, the parabolic spectroscopic response at each point over the surface. The other terms describe the dynamic ∆f variations occurring when scanning the surface. Considering the case of spectroscopic curves, the maximum of the ∆f vs. V DC curve is reached when:

V (c) DC = ∂∆f ∂V DC V mod =0 = 0 → V (c) DC = V cpd - C 1 2C 2 a 1 b 1 e +αz min (12)
Hence, the value of the compensated CPD in FM-KPFM is almost similar to the one measured in AM-KPFM. The difference stems from the ratios a 0 /b 0 and a 1 /b 1 . Expressions 10 and 12 show that the compensated CPD in AM-KPFM is sensitive to the averaged value of the electrostatic force (0-order Fourier component of the force), whereas it is sensitive to its Fourier component at the vibration frequency of the fundamental eigenmode in FM-KPFM.

To some extend, this illustrates the long-standing idea that AM-KPFM is sensitive to the electrostatic force and that FM-KPFM rather to its gradient 32 . Coefficients a n and b n only depend on the vibration amplitude of the fundamental eigenmode, A 0 . However, for n = 0 and 1, they weakly depend on it as soon as A 0 is larger than a few nanometers. With A 0 = 4 nm, one has a 0 = 0.1155, a 1 = 0.1080 and b 0 = 0.0724, b 1 = 0.0766. Hence FM-and AM-KPFM will not provide similar values for the local CPD, but they will neither differ significantly.

Numeric approach

The numeric implementation of the FM-KPFM setup is easier than the AM-KPFM one as the modulated signal is provided by the PLL without further processing. However, the sequence of the calculation of the spectroscopy curve slightly differs from the AM-KPFM one. The first part of the calculation is similar. When the approach is complete, we set V mod = 0, consistently with what is done during the experiments. Then, the DC part of the bias is swept. Similar sets of parameters have been used.

Results and discussion

Analytic and numeric ∆f vs. V DC spectroscopic curves are reported in fig. Finally, we would like to point out the following issue. Although the analytic variations of ∆f exhibit a quadratic dependence with V DC , readily visible in equ.11, it is impossible to fit the numerical variations with a quadratic polynom of V DC to a good accuracy (cf. fig. 4). So far, this behavior remains not understood. Although numerous attempts, including careful investigations of the influence of the PLL, we did not manage to identify the origin of this effect. It is obviously possible to fit the ∆f variations with a better accuracy when using higher-order polynoms, but this does not answer the question of the origin of such power laws, which do not stem from the expression of the electrostatic force. The above elements have demonstrated that the compensated local CPD was strongly sensitive to the tip's apex geometry (radii R and R a ). In order to assess how strong this dependence is, we have used equ.12 (FM-KPFM case) to plot the variations of the local CPD vs. R and R a at constant height (cf. fig. 5a). This is only made possible analytically, although the analytic method is less accurate than the numeric one. The parameters are I. For these tips, quadratic fits have been performed (continuous red and blue curves, respectively), which do not fit the numerical variations to a good accuracy (cf. text).

identical to those given before, in particular z min = 3 Å on top of an anion. The mesoscopic radius R is varied from 3 to 15 nm and R a from 1 to 1.6 Å. It is not allowed to make R a bigger owing to the main assumption of the model R a ≪ R. Despite the rather narrow range of R a , the local CPD is strongly decreased from -1.3 V (cf. also table I) down to -9 V as R a gets bigger. It is not possible to provide a tractable expression of the local CPD with R a owing to its complex and implicit dependencies (coefficients D 1 and A 2 ). However, some keys for understanding this behavior have been given in section IV A. On the opposite, the CPD does not exhibit a noticeable dependence with R, meaning that the nano-asperity plays a key role in the occurrence of the high-resolution imaging in KPFM. In order to extend the present framework, the influence of a long-range capacitive force on the measured local CPD has been assessed. We stick to the geometry defined before:

the counter-electrode, the bulk dielectric and the tip on top of it, but we consider an additional macroscopic metallic electrode placed a few microns above the tip [START_REF]Usual silicon-based commercial cantilevers have integrated tips with 10 microns height[END_REF] , hence several millimeters away from the counter-electrode. The electrode itself has millimetric dimensions and is at the same potential than the tip, V b . It depicts, to some extend, the cantilever and the tip holder when mounted in the microscope. With such a configuration, since the MSP decays exponentially with the tip-surface separation, it does not influence this electrode. Therefore, on the electrostatic level, the "mesoscopic/microscopic" problem and the "macroscopic" problem are rather decoupled, which is equivalent to include an additional long-range component to the total electrostatic force. Hence, equ.4 is modified into:

F es (V b , z) = F (1) mµ + F (2) µ + F M , (13)
To derive an approximate expression of F M , let's consider (1)-a planar capacitor with surface S and a distance between electrodes z M corresponding to several millimeters and (2)-that the dielectric sample fills entirely the vacuum between the electrodes. The condition 1 actually implies to neglect z, which is of the order of a few angströms, with respect to z M each time it is necessary. Then, we get:

F M (V b , z) ≃ - C ′ 2 (z M + z) 2 (V b (t) -V cpd ) 2 , ( 14
) with C ′ 2 = ǫ 0 ǫ 2 d S/2. A quick estimate for S = 1 mm 2 , z M = 5 mm (C ′ 2 = 8.5 10 -5 pN.m 2 .V -2 ) and V b = 1 V yields F M ≃ -3 pN.
For the FM-KPFM case, when applying Giessibl's formalism, an analytical expression for the frequency shift connected to F M can be derived. When applying a similar procedure than the one used to derive equ.12 to the frequency shift resulting from the total force (equ.13), the expression of the compensated CPD now becomes:

V (c) DC = V cpd - C 1 a 1 e -αz min 2C 2 b 1 e -2αz min + 4πC ′ 2 A 0 /z 3 M (15)
It is reminded that A 0 is the vibration amplitude of the first cantilever eigenmode. The graph of V (c) DC vs. z min is reported in fig. 5b for tips 1, 2 and 3. The diverging behavior predicted by equs.10 and 12 has disappeared because the term carrying the exponential decay in the denominator (deriving from F

(2) µ ) is now compensated by the long-range one. Therefore, the decaying exponential in the numerator (deriving from F

(1) mµ ) nullifies the local CPD far from the surface. This behavior is more likely to be measured and intricately originates from long-range electrostatic, macroscopic, interactions. The shape of the curve looks like the resonance curve of a harmonic oscillator. Hence, there is an optimum tip-surface separation for which, the strength of the KPFM contrast will be enhanced, although the measured value will not be correlated to the MSP. Upon the tip geometry, the maximum of the curve occurs at various tip-surface separations. In particular, the discrepancies between tips 1 and 2 suggest that it is more likely to achieve a stabler KPFM imaging regime with tips including a not too big mesoscopic apex. Let us notice that such a behavior is expected to stand also for the AM-KPFM case owing to the analogy between equs.10 and 12.

At last, we are aware that the predicted variations of the compensated local CPD are very large. This is the consequence of the power laws occurring in the expression of the short-range electrostatic force, themselves deriving from the raw geometry of the tip. That large variations will never be measured on the experimental level owing to realistic tip shapes which differ significantly from the one of our model. Nevertheless, we still believe that the reported variations are qualitatively relevant.

VI. CONCLUSION

The influence of short-range electrostatic forces on the measured local CPD by means of AM-and FM-KPFM methods has been discussed. The sample is a bulk alkali halide single crystal. This work relies on numeric and analytic descriptions of both KPFM modes. The analytic expression for the short-range electrostatic force is derived from a previous work and exhibits a short-range distance dependance and a linear as well as quadratic dependence of the bias voltage. The linear term is responsible for the atomic-scale CPD contrast, while the quadratic one accounts for the detected signal by the KPFM methods. The latter term stems from the dynamic polarisation of the ions of the crystal. Analytic and numeric approaches stress the influence of the linear term on the measured compensated CPD which makes it deviating not only from the actual value of the CPD, but also from the surface potential.

Hence, in the short-range regime, the AM-or FM-KPFM compensated local CPD neither reflects the tip-surface CPD nor the local surface potential, but rather an effective value which is convoluted by the geometry of the tip. Nevertheless, it is shown that the lateral periodicity of the signal is preserved, which makes the AM-or FM-KPFM methods indeed sensitive to atomic scale changes of the surface potential, as demonstrated experimentally, although quantitative numbers are unlikely derivable.

It has also been found that the local CPD detected by means of AM-or FM-KPFM is expected to differ. This stems from the detection methods which are not performed at similar frequencies, hence selecting different Fourier components of the modulated force for the signal processing. However, when using large amplitudes like those regularly used in beam-deflection based nc-AFM (a few nm), this effect is expected to be weak.

Finally, the influence of long-range, capacitive, electrostatic forces has been addressed in conjunction with the short-range ones. The local CPD then exhibits kind of a resonance peak as a function of the tip-surface separation. The position of the peak depends on the tip geometry, but is expected to occur at a few angströms above the surface. The occurrence of a resonance phenomenon illustrates the idea that there is an optimum tip-surface separation for achieving a large vertical KPFM contrast, although this contrast is not correlated to the surface potential.

3 .

 3 The set of parameters used to perform the calculations are consistent with the experimental conditions reported in19 : fundamental eigenmode: f 0 = 150 kHz, k 0 = 30 N/m and Q 0 = 30000; second eigenmode: f 1 = 936 kHz,k 1 = 1168 N/m, Q 1 = 8000.The modulation amplitude of the bias is V mod = 3 V. The oscillation amplitude of the fundamental eigenmode is A 0 = 4 nm (8 nm peak to peak). The curves have been computed for the three tips defined in fig.2aat z min = 3 Å , on top of an Comparison between analytic and numeric values of the local CPD (in volts) for tips 1 and 3 upon the KPFM method. The values have been estimated for z min = 3 Å above an anion of the KBr(001) facet. For this separation, the MSP is -440 mV. The CPD of the tip-surface system has been arbitrarily set to V cpd = 0 V.

FIG. 3 :

 3 FIG. 3: Analytic and numeric A 1 vs. V DC spectroscopic curves computed in AM-KPFM mode at z min = 3 Å with the three tips defined in fig.2. The values of the local CPD for tips 1 and 3 are given in table I.

V

  . OVERALL GEOMETRICAL INFLUENCE OF THE TIP TO THE LOCAL CPD: CONNEXION WITH LONG-RANGE ELECTROSTATIC FORCES A. Influence of the tip apex structure

FIG. 4 :

 4 FIG. 4: Analytic and numeric ∆f vs. V DC spectroscopic curves computed in FM-KPFM mode at z min = 3 Å with the three tips defined in fig.2. The values of the local CPD for tips 1 and 3 are given in tableI. For these tips, quadratic fits have been performed (continuous red and blue

B.

  Influence of long-range capacitive forcesSo far, we have been concerned about the influence of the electrostatic force stemming from the mesoscopic/microscopic structure of the tip on the measured local CPD, without taking into account long-range capacitive forces. Nevertheless, on the experimental level, these forces are always present. Furthermore equs.10 and 12 exhibit a surprising exponential dependence of the local CPD with the tip-surface separation, which should lead to a divergence/saturation of the measurement when attempting to perform it far from the surface.

FIG. 5

 5 FIG. 5: a-Evolution of the compensated local CPD as a function of tip's radii R and R a computed from equ.12, i.e. for the FM-KPFM case. A strong R a dependence is noticed. b-Compensated local CPD vs. tip-surface separation computed from equ.15 for tips 1 (black), 2 (dark grey) and 3 (light grey). The diverging behavior initially predicted when only considering short-range electrostatic forces is now disappeared. Instead, a resonant behavior is observed.

  4. For that set of simulations, we have used f mod = 2.5 kHz. The values of the local CPDs are reported in table I. As stated in the former section, the analytic values of the CPD slightly differ from those deduced by AM-KPFM. Similar trends are observed for the numerical results. The qualitative behavior between analytic and numeric curves is also similar, the measured local CPD with tip 2 being much shifted compared to the one measured with tips 1 and 3. Again, the local CPD's neither reflect the CPD of the tip-surface system, nor the MSP at similar separation.

Appendix

Out of ref. 19, the coefficients l 0 , l 1 , l 2 and l 3 can be written:

)