Jean-Claude Boulet

Jean-Michel Roger

Improvement of Direct Calibration in spectroscopy

Several linear calibration methods have been proposed for predicting the concentration of a particular compound from a spectrum. Some methods are based on experimental data, such as Partial Least Square Regression. Other methods are based on expert data, e.g. Direct Calibration. This article proposes a new method, called Improved Direct Calibration, which uses expert and experimental information. It performs a projection onto the pure interest spectrum, after correcting it from influence factors. No calibration dataset is necessary to build this model. This method has been successfully applied to the quantification of ethanol in musts during fermentation, using near infra-red spectrometry.

Introduction

Tests applied to many agrifood, pharmaceutical or chemical processes involve quantifying a factor of interest: for example, alcoholic fermentation is tested by measuring the ethanol produced by yeasts. Infrared spectroscopy tools are increasingly used for these applications. However, a spectrum does not only provide information specific to the factor of interest, it also contains contributions from influence factors such as the concentration of other chemical elements, temperature or turbidity. Nevertheless the value of the factor of interest can be extracted using chemometric tools. These build a calibration based on expert and/or experimental knowledge.

Expert knowledge is fundamental data regardless of the experiment. A pure spectrum, a chemical composition or a molecular weight are all examples of expert knowledge linked to a chemical factor. Some calibration methods use pure spectra, for example Direct Calibration (DC) [START_REF] Martens | Multivariate Calibration[END_REF], to produce a direct model. Extended Multiplicative Scattering Correction (EMSC) [START_REF] Martens | Light scattering and light absorbance separated by extended multiplicative signal correction, application to near infra-red transmission analysis of powder mixtures[END_REF] uses baseline models to correct for the effect of a physical influence factor, the scattering.

Experimental knowledge is represented by data collected from samples: spectra and associated quantitative and qualitative values. A calibration database, consisting of a set of spectra associated with the values of the 1 factor of interest is thus an example of experimental knowledge. A matrix of spectra acquired according to an experimental design in which the influence factors vary is another example. The most conventional regression methods such as Principal Component Regression (PCR) [START_REF] Martens | Multivariate Calibration[END_REF] and Partial Least Square Regression (PLSR) [START_REF] Wold | The collinearity problem in linear regression, the partial least square (pls) approch to generalized inverses[END_REF] use the calibration database to produce a forward-inverse model. Other methods use experimental knowledge complementary to the calibration database to correct for the effect of influence factors. For example, Independent Interference Reduction (IIR) [START_REF] Hansen | Pre-processing method minimizing the need for reference analyses[END_REF] uses spectra for which the factor of interest has zero value, to identify and remove useless space. External Parameter Orthogonalisation [START_REF] Roger | Epo-pls external parameter orthogonalisation of pls, application to temperature-independant measurement of sugar contents in fruits[END_REF] performs the same correction but with spectra acquired according to an experimental design in which only the influence factors vary.

The most commonly used chemometric calibration methods are based only on experimental knowledge. They are constrained by the calibration database management, which must represent the expected variability in future samples, for both the factor of interest and the influence factors. Thus, a robust model requires an often tedious and costly experimental phase. However, using a mixture of experimental and expert data should reduce this constraint. For instance, Science-Based Calibration (SBC) [START_REF] Marbach | A new method for multivariate calibration[END_REF] produces a projection using the pure spectrum of the factor of interest, weighted by the effects of influence factors identified in an experimental database. However, as proposed by Marbach, this method does not specifically take into account the pure spectra of chemical influence factors. The present article proposes a new method, mixing experimental and expert knowledge while taking these two types of informations into account. The first part provides a theoretical description of the method. The second part presents an application of monitoring alcohol fermentation using infrared spectrometry, and finally the third part presents and discusses the results.

Theory

Generally speaking, matrices are in bold upper case, vectors in bold lower case and scalars in normal characters. Vectors are listed in columns. The layout used for matrix elements depends on their nature: elements which are the same length as the spectra are in lines, those which are the same length as the number of individuals are in columns.

Let x(P, 1) be a spectrum acquired for a sample and y the value of an interest factor associated with this sample; k(P, 1) the pure spectrum of the factor of interest; K the pure spectra matrix (Q, P ) of the Q constituents other than the factor of interest and t χ , of dimensions (Q, 1), their concentrations. Let X G be a matrix of spectra acquired while the influence factors vary and the interest factor does not. Three approaches have been proposed to characterize X G . Hansen [START_REF] Hansen | Pre-processing method minimizing the need for reference analyses[END_REF] uses a set of samples in which the factor of interest is naught. Marbach [START_REF] Marbach | A new method for multivariate calibration[END_REF] uses a set of samples measured for the same value of the factor of interest, then centered. Roger [START_REF] Roger | Epo-pls external parameter orthogonalisation of pls, application to temperature-independant measurement of sugar contents in fruits[END_REF] uses a set of samples measured then centered at different levels of the factor of interest. Centering in the last two cases removes the spectral influence of the interest factor.

The factor of interest contribution to x is yk. Chemical influences are the spectral results of concentrations of all the chemical compounds other than the compound of interest in the sample analysed. If all the chemical compounds obey the Beer-Lambert law, each contributes a profile of its pure spectrum, weighted by its concentration, to the final spectrum. Their contribution to the spectrum is therefore given by: K t χ . The physical influence factors include all the disturbances involved in measuring the spectrum, such as temperature and granulometry for example. Even if they cannot be represented by pure spectra, their influences evolve in a subspace of R P . Let {p 1 ...p A } be an orthogonal basis in this subspace and P the matrix (A, P ) containing them. Let t φ be the score vector representing the physical influences on this basis. The contribution to the spectrum of physical values is therefore given by P t φ . x can therefore be written:

x = yk + K t χ + P t φ + ε
where ε is a noise vector, each component of which is independent of the others and identically distributed. Assuming that this noise is weak enough to be negligible, the previous expression leads to:

x = yk + K t χ + P t φ (1)
Calibrating a linear model consists in determining the b-coefficient vector b and the intercept b 0 such that y defined by:

y = x b + b 0 (2)
is the best estimate of y that minimises |y -y| under certain constraints. Several direct calibration strategies have been proposed, which give different results in the light of equation [START_REF] Martens | Multivariate Calibration[END_REF]. They differ in the way in which influence factor information is handled.

DC proposes to project x onto k orthogonally to K. Two conditions are assumed to be fulfilled: (c1) the pure spectra of all the chemical values are known and linearly independent so that (KK ) is invertible; (c2) the effect on the spectra of physical factors is assumed to be negligible. Let Σ DC of dimension (P, P ) be the K-orthogonal projector:

Σ DC = (I -K (KK ) -1 K) ( 3 )
Transposing and right multiplying by Σ DC the equation (1) yields:

x Σ DC = yk Σ DC + t χ KΣ DC + t φ PΣ DC
The two terms on the right of the equation are null, the first by construction, the second in application of the previous hypothesis (c2), and finally:

x Σ DC = yk Σ DC Right multiplying by k(k Σ DC k) -1 produces the DC formula if a single factor of interest is predicted [START_REF] Marbach | A new method for multivariate calibration[END_REF]:

y = x Σ DC k(k Σ DC k) -1 (4)
leading to:

b DC = Σ DC k(k Σ DC k) -1 (5)
Very attractive in principle, this approach is very difficult to apply because hypotheses (c1) and (c2) are rarely fulfilled.

SBC proposes to project x onto k using noise-reduction metrics. The SBC b-coefficients are given by:

b SBC = Σ SBC k(k Σ SBC k) -1 (6)
where Σ SBC is given by:

Σ SBC = (X G X G ) -1 (7)
Each variable of x is weighted. The weight of each variable depends on its variability on X G plus those of other correlated variables. Highest variabilities in X G lead to the lowest scores. This approach gives relatively more weight to noise-free variables of x which contain information about the interest factor. Therefore, SBC increases the signal-to-noise ratio. Problems can occur when X G is not invertible. Pseudo-inverse calculation or matrix dimension reduction using PCA are alternative solutions.

The new method presented in this paper (IDC) is performed by completing the expert information in K with the experimental information in X G . For this, a vector basis {p 1 ...p A }, representing the space spanned by X G , forming matrix P, is identified and added to K, giving a matrix R. Let Σ IDC be the R orthogonal projector:

Σ IDC = (I -R (RR ) -1 R) ( 8 )
By transposing, then right multiplying by Σ IDC the equation 1, the effect of physical and chemical influence factors becomes mathematically null, which gives:

x Σ IDC = yk Σ IDC Right multiplying by k(k Σ IDC k) -1 provides: y = x Σ IDC k(k Σ IDC k) -1 (9)
which gives:

b IDC = Σ IDC k(k Σ IDC k ) -1
The easiest way of identifying P is to perform an SVD on X G , and to retain the first A eigenvectors. The choice of A is an important stage in this method. On the one hand, if A is too high, residual information on the interest factor that could be present in X G can be captured by P, leading to a poor prediction. On the other hand, if A is too small, not all the influence factors are corrected, leading also to a poor prediction. To choose A, calculating the prediction error as a function of A onto a prediction set is possible. However, this approach loses the major advantage of direct calibration, which is precisely not to need a calibration database. Another approach consists of applying the IDC model to X G (for which the interest factor is constant and usually known) for different values of A and examining the evolution of the prediction error.

Although matrix K must theoretically contain the pure spectra of all the compounds in the sample, in practice this approach comes up against a certain number of difficulties. For instance, in very complex samples, some pure products cannot be extracted and stabilised in large enough quantities for spectra measurement to be possible. In order to capture all the information concerning influence factors in either X G or K, a simple rule should consist in considering that K must contain the pure spectra of products whose concentration does not vary in X G , or that the influence factors not found in K must be varied in X G .

Materials and methods

The application concerned ethanol quantification in musts and wines using near infrared spectroscopy. Spectra were acquired for all samples analyzed in the Skalli laboratory (Sète, France) during Septembre 2005, i.e. 80p.cent were fermenting musts, and 20p.cent were wines from previous vintages and sometimes other wineries, controled during the buying, ageing or bottling processes. Due to changes in the chemical composition of musts during vinification, all samples will be considered as different. The experimental database consisted of 1480 spectra, acquired using a Jasco spectrophotometer (optical path 1 mm, range 500 -2500 nm, step 2nm, with water as optical reference). The ethanol reference values of these same samples were measured by mid-infrared spectrophotometry (Foss). These data yielded an X matrix (1480, 1001) and a y vector (1480, 1). The pure spectra were acquired using the same spectrophotometer onto pure samples of ethanol, glycerol, water and lactic acid, with air as optical reference.

Data processing

The data were processed using Scilab software (www.scilab.org) completed with PCA and PLSR functions using the Saisir toolbox (www.chimiometrie.fr). All spectra baselines were previously adjusted to value 0 for wavelength 1100nm, chosen in a region where experimental and expert spectra both had the lowest variabilities. The interest factor was ethanol, its pure spectrum k was divided by 100 for the predicted values to be directly expressed in percent volume. Experimental information was represented by (X, y) and split into three sets:

• X 1 , containing non fermented musts, i.e. 165 sample spectra with zero ethanol concentration;

• (X 2 , y 2 ), containing the first 315 samples with non-zero ethanol content, in their chronological order of acquisition;

• (X 3 , y 3 ) containing the last 1000 samples with non-zero ethanol content, in their chronological order of acquisition.

This chronological splitting was chosen to ensure maximum independence between the sets of (X 2 , y 2 ) and (X 3 , y 3 ). It was verified that the histograms of y 2 and y 3 were comparable. The P matrix was obtained from the A first eigenvectors of an SVD of X 1 taken here as the X G matrix. Expert information about influence factors was represented by glycerol, lactic acid and water spectra, yielding to the K matrix (1001, 3). Seven calibration models were calculated then tested on (X 3 , y 3 ). The first three models were designed to explain how IDC works. The first model (m1) was a simple projection onto k. The second model (m2) used IDC with only k and K, corresponding to DC with very few pure spectra. The third model (m3) used IDC with only k and P. The fourth model (m4) used complete IDC with k, K and P. The fifth model (m5) used PLSR, calibrated on (X 2 , y 2 ) by cross-validation of the NIPALS algorithm. The number of latent variables was chosen to minimise the RMSECV. The sixth model (m6) was an IDC with k, K and P with A distinctly higher than the optimal value chosen in models (m3) and (m4). Finally, the seventh model 6

Author-produced version of the article published in Analytica Chimica Acta, 2010, 668, 2, 130-136. The original publication is available at http://www.sciencedirect.com/ doi:10.1016/j.aca.2010.04.010 (m7) repeated model (m4) after the water spectrum had been removed from K.

Model comparison

The models were first evaluated visually according to their general predictive aptitude, i.e. by aligning the predicted values relative to the reference values along the line ( y = y). RMSEP, bias and RMSEPc corrected for bias were calculated for each model. The b-coefficient peaks were interpreted by comparison with the pure spectra of ethanol and water.

Results

Choice of y unit

The Beer-Lambert law is limited to low concentrations, which is not the case here. However, absorbance is due to electron excitement by photons, so y should reflect electronic properties of the molecules. Mark [START_REF] Mark | Chemometric calibration without matrices (almost)[END_REF] recently raised the question of the unit of predicted values using DC. He suggested that y represent the percentage of hydrogen atoms contributed by the compound of interest. The number of H-moles contributed by 1mL of ethanol and water is almost the same, respectively 0.104 and 0.111, so in first approximation y would represent the percent volume of ethanol, that is the common unit used in enology.

The unit of y is given by the unit chosen for k. Here, k was the pure spectrum of ethanol, so k/100 represents the spectrum of 1 p.cent vol. of ethanol, and the corresponding y values obtained using the IDC model are also expressed as percent volume of ethanol.

Construction of the K matrix

It is unusual, in spectroscopy, to use spectra not acquired under the same conditions. However, this is the case here, because the experimental data represented by X 1 , X 2 and X 3 and the expert data represented by k and K used water and air respectively as optical references. The explanation is that experimental spectra are usually acquired with a water reference for practical reasons, whereas expert data are acquired with an air reference to meet the requirements of Beer-Lambert's law relative to mixtures. For any spectrum, let note x w its value with a water reference and x a its value with an air reference. If k a water is the water spectrum with an air reference, then: Therefore, the difference between the different x spectra measured with air and water references is simply the water spectra. Once this same spectrum is inserted in K, because Σ IDC is an orthogonal projection onto K, the product Σ IDC k a water is null. In conclusion, incorporating the water spectrum k a water in K means that the experimental data acquired with water reference can be used directly instead of the experimental data with air reference.

x w = x a -
The main natural constituents of musts and wines, excluding ethanol, are: water, glucose, fructose, glycerol, tartaric, malic and lactic acids. However, the musts used for X G , named X 1 in our application, contain variable quantities of glucose and fructose as well as tartaric and malic acids. The pure spectra of these compounds were therefore not included in K. On the other hand, glycerol and lactic acid are not found in musts, so they were not present in X G . This is why their pure spectra were included in K. Finally, K contains spectra for water, glycerol and lactic acid.

Parameter determination

Fig. 1 is used to choose the dimensions of models (m3), (m4), (m6) and (m7). Fig. 1a shows the evolution of the percentage of inertia of X G captured by P vectors. Fig. 1b shows the evolution of the standard error of prediction for model (m4) applied to X G as a function of A. The value A = 4 was used because it allowed almost all X G information to be captured while giving a minimal prediction error. This optimal value of A = 4 was also applied to (m2) and (m7). The prediction error increases for A > 10. That confirms the risk that excessively high values of A can cause P to capture residual information about the factor of interest, as mentioned in the theory part. To verify this, model (m6) with A = 12 was also constructed. Fig. 2 was used to tune the dimension of the PLS model. The standard cross-validation error of the PLSR (RMSECV) was stabilised for 5 and minimum for 8 latent variables, so the PLSR model was built using 8 latent variables.

Analysis of models m1 to m7

Fig. 3-m2 shows the prediction obtained by model (m2). The correction induced by the pure spectra led to a prediction which certainly includes too much noise to be used as a prediction model, but which is sensitive to y: the square correlation coefficient between predicted values and reference values is 0.87.

Table 1 shows the correlations obtained by different models, derived from (m2) by eliminating 0 to 2 spectra from the K matrix. Highest correlations are clearly obtained by the presence of the pure spectrum of water and at least one of the other two pure spectra: glycerol or lactic acid. In all other cases, correlation does not exceed 0.20 and can be considered to be null. The need of the water spectrum is explained by the different references used for data acquitition, air or water. The cosines of glycerol and lactic acid spectra with the ethanol spectrum are 0.93 and 0.92 respectively. This strong colinearity removes from the spectra a huge amount of useless information close to the ethanol signal. The remaining information is more specific to ethanol, signal to noise ratio is increased, leading to more accurate models.

Fig. 3-m3 shows the prediction obtained by the model (m3). This model does not use any of the three spectra mentioned previously: water, lactic acid, glycerol. It is therefore normal for the correlation between predicted and reference values to be close to 0. However, the introduction of X G eliminates all prediction variability and thus noise. Finally, this model predicts a value close to 0 for all the samples. Fig. 3-m1 shows the prediction obtained with no correction. The absence of K leads to a prediction which is insensitive to ethanol, and the absence of P gives a very noisy prediction. Logically, the prediction for (m1) is very noisy and centered on 0.

The complete IDC model (m4) gives very satisfactory predictions, comparable with those of the PLSR (m5), see Fig. 3-m4 and m5, Tables 2 and3. The IDC RMSEP is very close to that of the PLSR for ethanol concentrations of less than 10p.cent vol., but higher for ethanol concentrations over 10p.cent vol. In general, models (m2) to (m4) display greater error in the high ethanol concentration zone, corresponding to wines at the end of fermentation or finished wines. The problem is certainly due to the appearance of new influence factors, such as an effect of physical and chemical stabilisation of finished wines. The PLSR was able to handle these new influence factors because the calibration set contained finished wines, thus improving robustness for PLSR in this situation. These new influence factors haven't been taken into account by IDC calibration because they were represented neither in K nor in X G . Another set of data from finished wines could have been be taken to improve IDC calibration. The problem is then to obtain samples which all have exactly the same ethanol content so that simply centering eliminates the effect of ethanol on these spectra. A solution can be an experimental design based on samples before and after stabilisation, under such conditions that no loses of ethanol are possible.

Model (m6), Fig. 3-m6 had poor prediction, much worse than that of model (m4). The difference between these models is that the optimal value A = 4 was chosen for (m4), whereas a value deliberately chosen to be much higher, A = 12, was chosen for (m6). Thus, the experimental design leading to X G was intended to maximize the expression of the influence factors, which are modelised by the first loadings of P. But native grapes can also contain low amounts of ethanol, because of anaerobic metabolism known as carbonic maceration, or yeast activity on wounded berries. This ethanol can lead to a weak spectral signal. If the dimension of P is too large, more or less information about the interest factor is captured by X G . Thus, Σk tends to the null vector. This was verified in our example, an IDC model with A = 40 led to predictions centered around 0, results not shown.

The difference between models (m4) and (m7) lies only in the presence of the water spectrum for the first one, and its absence for the second one. Their differences of accuracy demonstrates the importance of the water spectrum into K.

Fig. 4 represents the IDC b-coefficients compared to pure spectra of ethanol and water, completed by a wine spectrum. The 4 main peaks of the ethanol spectrum (1580, 1710, 2085 and 2295 nm) are found in the IDC b-coefficients (Fig. 4). The ethanol peak at 2085 nm is attenuated in the b-coefficients. One explanation is that sugars also have a strong absorbance peak at this wavelength [START_REF] Osborne | Near infrared spectroscopy in food analysis[END_REF]. In addition to the ethanol peaks, two other peaks are visible in the b-coefficients. The negative peak at 1450 nm as well as the positive peak at 1940 nm may be linked to the water absorbance in this zone. For ethanol prediction using IDC, the opposite contributions of these two peaks cancelled out the influence of water, that is precisely the interest of the orthogonal projection.

The comparison between PLSR b-coefficients and IDC shows that they are different: their cosine is 0.47, so the angle between these two vectors is nearly 60 degrees. This example illustrates the non-uniqueness of the models: equivalent predictions are obtained from very different models.

Bias and slope management

From a theoretical point of view, a well-built IDC model do not have bias and slope different from 0 and 1 respectively. Otherwise that means that some influence factors haven't been taken into account.

An IDC model can be built in certain conditions then applied under conditions where an unexpected influence factor hasn't been taken into account Let us assume that this new influence factor is described by only one vector e. For a sample i, let x i and y i be respectively the spectrum and the IDC estimation that would be obtained without the new influence factor, x * i and y * i the measured spectrum and the IDC estimation obtained with the new influence factor. Then there exists a scalar λ i that verifies:

x * i = x i + λ i e (10)
and after multiplying each side by b IDC :

y * i = y i + λ i b IDC e (11)
This explains why some IDC models present a non-null biais, a slope different from 1, or both. When the {λ i } coefficients take a constant value λ for all i, predicted values have a constant biais equal to λb IDC e. When the {λ i } coefficients are correlated to the y i , the resulting predictions have a slope different from 1. In the other cases, the prediction variance is increased.

The bias problem can be due to the way X G is built. A first approach consists in acquiring spectra with a constant value of y and variation of influence factors as proposed by Roger [START_REF] Roger | Epo-pls external parameter orthogonalisation of pls, application to temperature-independant measurement of sugar contents in fruits[END_REF] and Marbach [START_REF] Marbach | A new method for multivariate calibration[END_REF]. The resulting matrix X G must be centered to withdraw all information about the interest factor. This operation also eliminates all constant information, so that the correction of any constant influence factor e is not performed. A second approach proposed by Hansen [START_REF] Hansen | Pre-processing method minimizing the need for reference analyses[END_REF] uses a set of samples where the interest factor is null whereas influence factors vary, as used in this paper. No centering is necessary, so the effect of e is taken into account.

However if e is known, it is always better to put it in K in order to take it into account explicitely. This can be illustrated by a simple example. Let x G be a vector of X G . The orthogonal projection to the vector (x G +e) removes a space of dimension 1, whereas the orthogonal projection to the matrix containing x G and e removes a space of dimension 2. All the information about x G and e is withdrawn in the second case, but not in the first. In our application, e is the water spectrum used to get spectra with water reference. Despite its presence in all spectra of X G , it must be put into K to obtain an accurate IDC model.

Links with the Net Analyte Signal

The term (Σ DC k) is written (I -K (KK ) -1 K)k; it represents the projection of the pure spectrum of the factor of interest orthogonally to the matrix of pure spectra of chemical influence factors. Similarly, (Σ IDC k) is written (I -R (RR ) -1 R)k. Therefore the k spectrum of the factor of interest is projected orthogonally to the space defining the chemical and physical influence factors. These two cases are in accordance with the definition of the Net Analyte Signal (NAS) given by Lorber [START_REF] Lorber | Net analyte signal calculation in multivariate calibration[END_REF]:"the net analyte signal may be computed as the part of its spectrum orthogonal to the contribution of other coexisting constituents", with the difference that with IDC this definition is extended to physical influence factors: IDC improves the definition of NAS. Let the scalar α = (kΣ IDC k ) -1 . The prediction of the interest factor is also written: y IDC = αX NAS IDC Therefore, prediction by IDC is proportional to the inner product of X spectra and the estimate of NAS calculated by IDC using Euclidean metrics. Coefficient α adjusts the notation scale of the interest factor which is arbitrary: for example, mg/L or g/L. IDC b-coefficients tend towards the NAS, which is not the case for PLSR b-coefficients.

Conclusion

This study confirms the fact that DC cannot be applied in cases where certain chemical and physical influence values are not taken into account. However, these expert data can be judiciously completed with experimental information acquired in the form of a spectral matrix which is used to characterise the influence factors missing in the expert data. Simultaneous use of these expert and experimental data leads to the IDC. IDC does not require a calibration set, but results may need bias and/or slope correction before use. It is also shown that IDC is a predictive method based on the NAS, i.e. the predicted value is proportional to the inner product between the sample spectrum and the NAS.

IDC is much more efficient than DC. Sometimes it could perform as well as PLSR. It should be noted that IDC and PLSR b-coefficients are not similar, even when their respective predictions are close to each other. The applications of IDC are mainly those in which PLSR is not applicable or can only be used with difficulty, i.e. when construction of the calibration database is problematic or impossible. For instance IDC should be a powerful tool when applied to hyperspectral image analysis. In these situations, IDC provides high analytical added value.

Indirectly, this study showed that three methods of direct calibration: DC, SCB and IDC are built on the same formula and based on two matrices: 1) the pure spectrum k of the factor of interest, and 2) a Σ matrix characterising the influence factors, defined differently by DC, SBC or IDC. In a future article, we shall show that this same writing provides links between direct calibration and inverse calibration methods such as PLSR or PCR.

8Figure 1 :Figure 2 :

 12 Figure 1: (a) Evolution of the percentage of inertia of X G captured by A for the first 15 vectors of P; (b) Standard error of prediction for the IDC model applied to X G for A = 1 to 15

Figure 4 :

 4 Figure 4: From top to bottom, shapes of pure ethanol spectrum, ICD bcoefficients, pure water spectrum and example of wine spectrum

  Author-produced version of the article published in Analytica Chimica Acta, 2010, 668, 2, 130-136. The original publication is available at http://www.sciencedirect.com/ doi:10.1016/j.aca.2010.04.010

k a water 7

11

  Author-produced version of the article published in Analytica Chimica Acta, 2010, 668, 2, 130-136. The original publication is available at http://www.sciencedirect.com/ doi:10.1016/j.aca.2010.04.010

	Spectra in K R 2
	Water (W)	0.20
	Lactate (L) 0.00
	Glycerol (G) 0.06
	L + G	0.03
	W + L	0.74
	W + G	0.85
	m1 W + L + G 0.87	m2
		m3	m4
		m5	m6

m7 Figure 3: Model tests: (m1) simple projection without correction; (m2) IDC using only K; (m3) IDC using only X G ; (m4) complete IDC with A = 4; (m5) PLSR using 8 latent variables; (m6) complete IDC with A = 12; (m7) model (m4) after withdraw of water spectrum

Table 1 :

 1 Coefficients of correlation R 2 between predicted and reference values, for models obtained from (m2) after removing 0, 1 or 2 spectra

	Model	Slope	Bias RM SEP c RM SEP	R 2
	m1 (Σ = I)	-0.34 11.3	9.98	15.1	0.02
	m2 (Σ = K)	1.17	11.9	1.87	12.1	0.87
	m3 (Σ = P, A = 4)	-0.16 11.6	4.55	12.4	0.29
	m4-IDC (Σ = [KP], A = 4)	0.99	0.03	0.96	0.96	0.94
	m5-PLSR (8LV)	0.97	0.09	0.85	0.85	0.95
	m6 (Σ = [KP], A = 12)	0.43	6.32	2.41	6.76	0.74
	m7 (m4 without water spectrum) 1.10 -1.51	1.51	2.13	0.89

Table 2 :

 2 Figures of merit of the models

	Model	Ethanol < 10 % Ethanol ≥ 10 %
	IDC (m4)	0.87	1.01
	PLS (m5)	0.85	0.86
	Table 3: RMSEP of IDC and PLSR according to the concentration of ethanol
	in the sample
			13

Author-produced version of the article published in Analytica Chimica Acta, 2010, 668, 2, 130-136. The original publication is available at http://www.sciencedirect.com/ doi:10.1016/j.aca.2010.04.010

A Acknowledgements

This work was made possible within the IRVIN program, carried out with the Skalli winery thanks to a financial support from the Languedoc-Roussillon Region.

B References

Author-produced version of the article published in Analytica Chimica Acta, 2010, 668,


Franck Bocquet

Laurent Nony

Christian Loppacher
email: christian.loppacher@im2np.fr.

Polarization effects in Non-Contact Atomic Force Microscopy: a key to model the tip-sample interaction above charged adatoms

Keywords: numbers: 07.79.Lh, 41.20.Cv, 34.20.-b noncontact-AFM, atomic resolution, short-range forces, electronic polarization, charge identification

We discuss the influence of short-range electrostatic forces, so called dipolar forces, between the tip of an atomic force microscope (AFM) and a surface carrying charged adatoms. Dipolar forces are of microscopic character and have their origin in the polarizability of the foremost atoms on tip and surface. In most experiments performed by non-contact AFM, other forces such as binding forces dominate the interaction. However, in the experiments presented by Gross et al. [Science 324, 1428 (2009)], where the charge state of individual gold atoms adsorbed on a thin dielectric layer was determined, binding forces are negligible as the tip-sample distance is relatively large.

We develop a model which mimics the experimental tip-sample geometry of the aforementioned experiments. The model includes van der Waals and long-range electrostatic interactions, as well as the short-range electrostatic interaction based on the self-consistent description of electronic polarization effects on neutral and charged adatoms. The model is based on a calculation of the electrostatic energy of the tip-sample geometry.

Our calculations of non-contact AFM imaging as well as of bias spectroscopic curves are in good agreement with the experimental ones presented by Gross et al. It is demonstrated that the short-range dipolar force is mainly responsible for the contrast observed in topography imaging above charged species. However, it is the long-range capacitive force which is responsible for the detection of the charge state in bias spectroscopy. We discuss implications of our findings on future experiments which aim to detect single charges by means of Kelvin probe force microscopy.

I. INTRODUCTION

Atomic force microscopy (AFM) 1 in the so called non-contact mode (nc-AFM) 2,3 was used for the first time by Giessibl 4 to measure the semiconducting surface of Si(111)7x7 with atomic resolution. Since then, nc-AFM developed to an almost standard tool for imaging not only semiconducting, but also metallic and insulating surfaces with atomic resolution (cf. Ref. [5]). Although nc-AFM imaging with atomic resolution is now routinely achieved, there is still a debate about the microscopic origin of the contrast formation.

For semiconducting surfaces there is an agreement on the fact that the contrast is basically obtained due to the covalent bonding between the foremost tip atom and the surface atoms.

For semiconducting surfaces with well known adsorbates, the distance dependence of these covalent binding forces can, for example, be used in comparison with theoretical calculations to chemically identify the adsorbed atom species [START_REF] Sugimoto | [END_REF] . In other experiments with very small tipsample separation, the angular component of the binding forces contributes to the contrast formation 7,8 , which results in subatomic features in the images. However, also for larger tip-sample separations, atomic scale contrast has been observed on Si(111)7x7 Ref. [9]. In this work, the authors use a non-reactive oxidized silicon tip and explain that the tip-sample interaction is dominated by the dipole moment induced in the tip due to the surface charge distribution and not by covalent bonding.

For metal surfaces, the origin of the atomic contrast in nc-AFM imaging is less clear. Some works point out that for both, a clean and thus reactive silicon tip as well as for a tip contaminated by a metal atom, the tip-sample interaction would be dominated by strong covalent binding. However, this can not explain several experimental results where only weak orbital overlap was suggested and the decay length of the force was two times larger above the surface than for the tunneling current 10,11 . These results could only be modelled by using a tip carrying a water molecule, which resulted in only weak interaction due to van der Waals (vdW) forces, dipole-dipole interactions, and charge transfer 12 .

In the case of insulating surfaces, the interaction forces are multifold and strongly depend on the system under investigation. Most of the work done on insulators concentrate on two main issues: the polarity of the tip as well as the influence of dipolar interaction. It was Giessibl that first proposed an electrostatic imaging mechanism for ionic crystals based on a polarizable tip 13 . The proposed decay length of the interaction force in this model is about one order of magnitude per Å and thus sufficiently small to explain atomic scale contrast on ionic surfaces. In a series of papers using different model tips on ionic as well as oxide surfaces, Shluger et al. point out the importance of the polarization of both, tip and sample surface [14][15][16][17][18] . These extensive calculations consider, in particular, polarization effects; however, it seems difficult to separate them from the different other interactions also described in these works. The polarization is induced by the electric field of the tip.

The most prominent effect due to a polar tip is the local polarization of the surface which goes along with (vertical) displacements of surface atoms. The importance of knowing the polarity of the tip when interpreting the imaging mechanisms on ionic surfaces has been discussed in detail by different groups 19,20 .

When it comes to imaging of local charges and surface potentials, Kelvin probe force microscopy (KPFM) 21,22 is the method used in most experimental setups. KPFM is a method working in combination with nc-AFM and which is based on the minimization of the electrostatic interaction between the tip and the sample. Electrostatic tip-sample interactions have been widely studied in the long-range regime on the experimental and on the theoretical level [23][24][25][26][27][28][29] within the frame of nc-AFM and KPFM [30][31][32][33][34] . Long-range electrostatic forces cause capacitive forces connected to the capacitance of the tip/counter-electrode capacitor as well as Coulombic forces connected to the presence of charges within it.

For conducting surfaces it is well established, that KPFM reveals the so called local contact potential difference (LCPD) with a lateral resolution in the order of the tip front end [35][36][37][38][39][40] . For insulating surfaces and nanoscale contrast, the contrast mechanisms are much more complex since not only the exact geometry of the tip-sample system has to be considered, but also charges and (induced) dipoles within the system (polarization effects) [41][42][43] .

Electrostatic forces in the short-range regime are known to occur since atomicallyresolved LCPD images have already been reported on semi-conducting surfaces: Si(111)7×7

Refs. [31,44-47], Si(111)5 √ 3 × 5 √ 3-Sb Ref. [48], GaAs(110) Ref. [49] and InSb(001) Ref. [50],

and on two bulk dielectric surfaces: TiO 2 (110) Refs. [41,42] and KBr(001) Ref. [43]. However, the microscopic origin of these short-range electrostatic forces and the way they are coupled to the nc-AFM/KPFM imaging remains unclear, which has triggered intensive experimental and theoretical efforts in that field. Some groups initiated theoretical studies on the contribution of short-range electrostatic (SRE) forces 43,47,50,51 in nc-AFM/KPFM combined experiments. In a recent series of papers 43,52,53 forces may develop between a biased metallic tip carrying an atomic asperity such as a single atom or a small ionic cluster and the (001) facet of a bulk alkali halide single crystal. The results point out that these originate from the dynamic polarization of the atomic cluster and of surface atoms owing to the combined influence of the biased tip and of the Madelung surface potential of the ionic crystal 13 . In Ref. 57 where triboelectrification was measured with a sensitivity of a few charges only, single charges [58][59][60][61] , dipole moments [62][63][64] as well as the charge state of individual point defects 65 have been detected.

Due to their fundamental physical origin, SRE forces are ubiquitous. They are more and more under discussion since in recent years nc-AFM and KPFM are used to investigate single charges and local surface potentials on the nanometer length scale. In this work, we focus on the SRE interactions and ask two questions in regard to the imaging and the spectroscopy of charged species:

• can SRE forces alone be responsible for the atomic scale image contrast?

• is single charge detection on the nanometer scale possible?

These two questions have been raised by recent experimental results presented by Gross et al. 56 . The authors present a set of experimental results combining single electron charge sensitivity and atomic lateral resolution the latter, however, in topography imaging and not in KPFM. The investigated system is a charged or neutral gold atom adsorbed on a thin dielectric layer (NaCl) on a metal surface. The tip carries a foremost gold atom.

Different charge states of the adatoms are clearly resolved in nc-AFM topography images and bias spectroscopy can be used to identify the charge state of the gold adatom. Note that for the reported experimental tip-sample distances (> 4.5 Å), any interaction due to covalent binding can be excluded and only vdW and electrostatic forces have to be regarded.

Therefore, these experiments are an ideal model case in order to study the role of dipolar forces in nc-AFM experiments.

In the present work, we clarify the role of SRE interaction in the topography imaging mechanism of single charges and we give the limits for lateral resolution in single charge detection. Thereto, we develop a relatively simple analytical model which mimics the experimental tip-sample geometry as presented by Gross et al. and we use it to calculate the experimental images (line profiles) as well as the spectroscopy curves. In addition to the well known vdW forces [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF] and the long-range electrostatic forces 23,[START_REF] Guggisberg | [END_REF] for a known tip geometry, our model includes a self-consistent description of the SRE interaction including electronic polarization effects. Our purpose is to show how SRE contribution can reproduce trends and orders of magnitude reported in the experiments by Gross et al. 56 . However, it is not our intention to fit their results since that would require a fully atomistic description involving a complex simulation process [14][15][16][17][18] .

FIG. 1: Geometry of the problem. The tip consists of a macroscopic metallic body with a conical shape. The tip apex is described by a sphere with a radius R. The tip is grounded with regard to the counter-electrode on top of which lies a thin dielectric consisting of two mono-layers of an ionic film. The tip apex carries a neutral, but polarizable single atom. On top of the dielectric lies an adatom that may carry a single electron charge. The geometrical parameters of the problem are defined similarly to those in the work by Gross et al. 56 .

II. ELECTROSTATIC MODEL

A. Framework

The geometry of the problem is shown in Fig. 1. It accurately mimics the geometry of the problem as defined in the work by Gross et al. 56 . Our parameters are defined similarly as well. The tip consists of a metallic macroscopic body with a conical shape. The apical part of it is depicted by a sphere with a radius R. The apex carries a neutral single metallic atom in topmost position. The atom has an electronic polarizability α. The tip is grounded with regard to the metallic counter-electrode located below a thin dielectric with a thickness h d and biased at a potential V b . In Ref. [56], the dielectric is made of two mono-layers of NaCl deposited on a Cu(111) substrate. Thus, typically h d 3 Å. On top of the dielectric lies a metallic adatom that may carry a net single electron charge q = q × e, e being the elementary charge unit (1.6 × 10 -19 C). Hence, q = ±1, or q = 0. It is assumed that the surface adatom has the same electronic polarizability α as the one present on the tip apex.

Consistently with Ref. [56], the tip atom and the adatom are supposed to be Gold atoms. We therefore use the corresponding polarizability for that metal. In Ref. [68], pseudopotential calculations give an estimate for the dielectric susceptibility of Gold χ d = α 0 : 0.04×10 -39 < χ d (F.m 2 ) < 0.09 × 10 -39 . Thus, we will consider in the following χ d = 0.06 10 as in the present situation. Dipolar forces, as we define them, have a microscopic origin and are connected to the ability of the atoms at the interface to get polarized owing to the dual influence of the biased electrodes and of the presence of the charge within the tip/counterelectrode capacitor. Thus, the tip atom and the surface adatom carrying the charge q may acquire net dipolar moments, -→ p tip and -→ p ad , respectively. It is important to notice that these forces have a significantly different microscopic origin compared to vdW forces. This statement is justified in the following.

These forces have distinct decay lengthes, which permits to split their influence between the long-range (LR) regime (z > 1 nm) and the short-range (SR) regime (z < 1 nm).

Therefore, the total force acting between the tip and the surface may be written as:

-→ F tot = -→ F LR + -→ F SR = -→ F el LR + -→ F vdW LR + -→ F el SR + -→ F vdW SR (1)
The origins of LR-and SR-vdW forces are well-established. They are very briefly discussed in the following sections. At this stage, we focus on the electrostatic force, especially on the SR part of it, which has been under discussion in a few recent publications, as detailed in Sec. I.

In the present work, an expression of the SRE force regarding the system depicted in Fig. 1 is derived. The starting point of the approach is the expression of the electrostatic energy for the complete tip/vacuum/charge/dielectric layer/counter-electrode interface, W el described in Eq. ( 2). The vertical component of the total electrostatic force acting on the tip F el tot = F el LR + F el SR will then be derived as F el tot = -∂W el /∂z. The constitutive equation of the electrostatic energy for the system under investigation is

W el = - 1 2 C 0 V 2 b + qφ 0 + 1 2 qφ im LR + - 1 2 -→ p tip • -→ E ext tip - 1 2 -→ p ad • -→ E ext ad SR (2)
The first bracket stands for the usual electrostatic energy of the tip/vacuum/dielectric layer/counter-electrode capacitor including the inner charge, as described for instance in

Refs. [70,71]. In our problem, this term has a macroscopic origin since it depends on the geometry of the capacitor. Therefore, forces with a LR character will be associated with it. As said before, this term may be interpreted as the superposition between a purely capacitive force and Coulombic forces. C 0 depicts the capacitance of the former capacitor including neither the inner charge q nor the microscopic dipoles -→ p tip and -→ p ad . This defines a void capacitor, the capacitance of which may be derived upon integration over its geometry. φ 0 is the electrostatic potential which develops within the void capacitor at the location of the surface adatom, that is, the location of the charge q. φ im is the electrostatic potential induced by the set of image charges of q within both metallic electrodes at its location.

Although usual in electrostatics [START_REF] Jackson | Classical Electrodynamics[END_REF] , the second bracket in Eq. ( 2) is less usual when dealing with nc-AFM. It depicts the electrostatic energy of microscopic dielectrics, here the tip atom and the charged adatom, when polarized by an external field. In the present case, on the microscopic level, we may expect that the adatom gets polarized by the electric field which develops between the tip and the counter-electrode referred to as -→ E ext ad .

-→ E ext ad results from the superposition of two components: the field within the void capacitor plus the field induced by the image charges of the charged adatom within the metallic electrodes. Thus, the adatom will acquire a dipolar moment -→ p ad . Simultaneously, the tip atom gets polarized (dipolar moment -→ p tip ) by the electric field -→ E ext tip that may be built in the same manner: the field of the void capacitor, the field induced by the image charges of the charged adatom within the electrodes, but also the Coulombic field developed by the charge itself.

-→ E ext ad does not include the latter term since a charge can not influence itself.

Thus, the expressions of the microscopic dipoles -→ p ad and -→ p tip and those of the external electric fields

-→ E ext ad and

-→ E ext tip are to be derived. For that purpose, few assumptions are required.

B. Main assumptions

It is assumed that the dipoles -→ p ad and -→ p tip have a microscopic origin because they are only connected to the electronic polarizability of the adatom and of the tip atom, α. Hence, any depolarizing field is considered, as it should be the case with an extended dielectric.

In order to keep the model as simple as possible, it is retstated that the tip atom and the surface adatom consist of the same material. Therefore they have the same polarizability α, which is not strictly correct if the adatom is charged.

The influence of the thin dielectric layer within the capacitor is described by means of an effective dielectric permittivity ef f , instead of 0 . This would be correct if one was dealing with a planar capacitor, which is false due to the shape of the tip apex. However, we show that dipolar forces have a SR character. Therefore, the strongest effects are expected to occur when the adatom is placed right below the tip atom at a distance of a few Angströms.

Then, the "local capacitor" influencing the adatom is mainly localized within an almost planar area corresponding to the foremost part of the spherical apex [START_REF]A similar[END_REF] . This assumption is valid for R z. The capacitance of a planar capacitor, the electrodes of which are at a distance of h from each other, is C = 0 /h. When including a dielectric with a permittivity d and a thickness h d between the electrodes, the capacitance becomes:

1/C = (h -h d )/ 0 + h d / d . Therefore: C = ef f /h, with ef f = 0 d /[ d (1 -h d /h) + 0 h d /h].
Considering parameters corresponding to the geometry of the problem: h 8 Å, h d 3 Å, d = 5 0 (NaCl), ef f 1.4 0 . Thus, the influence of the discontinuous dielectric within the capacitor consists in replacing 0 by an effective dielectric permittivity with a slightly larger value. The preceding assumption implies that the dielectric is treated as a continuous medium. Hence, no microscopic description of it is considered.

Following the preceding elements, the polarizing fields influencing the adatom and tip atom referred to as -→ E dip ad and

-→ E dip tip , respectively (cf. Sec. II D), are calculated by assuming a planar geometry of the tip, which implies R z as well.

We do not address the influence of quantum effects in the problem and stick to a classical description. Hence, we attempt to describe neither the interaction between the adatom and the ionic layer nor the stability of the tip atom at the apex and the way its polarization might influence the local electronic density of the apex, thus modifying the resulting dipole, by screening effects for instance. In this approach, the interfacial atoms are somewhat frozen and the resulting interaction felt by the tip is measured at distances much larger than the tip atom/apex or adatom/ionic layer distances. These assumptions look reasonable in that the experimental results by Gross et al. are performed for tip-surface separations larger than 4 Å and at 5 K, which ensures the atomic stability of the tip apex and of the adatom over the surface. Nevertheless, from a classical point of view, the polarized tip atom creates an image dipole within the metallic apex, very close to the source dipole, with the same strength as it. Thus, because those dipoles are very close to each other, we assume that this ultimately produces an equivalent dipole built as their vectorial summation (cf. Appendix).

Because the problem is being treated in a classical manner, it is assumed that no wave function overlap between the electrons of the tip atom and those of the adatom takes place.

In other words, so called "chemical SR interactions" are not taken into account. This implies that the tip-surface separation remains larger than 4-5 Å, which is indeed the case in the work reported by Gross et al. At these distances, this also implies that the Madelung surface potential of the NaCl thin film is too weak to significantly modify the local electrostatic potential, owing to its exponential-like decay 43,[START_REF] Watson | [END_REF] . Thus, the influence of the Madelung potential on the resulting SRE force is not described. As a consequence, it must be noticed that the present description significantly differs from the models by Bocquet et al. 43 and Giessibl 13 on a bulk alkali halide crystal, for which polarization effects uniquely stem from the Madelung surface potential.

For the sake of clarity, the way the analytical expressions are derived is not detailed in the text. We rather discuss the main issues deriving from the formal expressions. The details of the analytical method and the required calculations are extensively detailed in the Appendix.

C. Long-range part of the total force

The preceding elements give the following expression for the force acting between the tip and the surface in the LR regime:

F LR = F vdW LR + F el LR = F vdW LR - ∂ ∂z - 1 2 C 0 V 2 b + qφ 0 + 1 2 qφ im (3)
F vdW LR depicts the London dispersion force between the tip and the surface. It consists of two contributions. The first one is calculated in a cone+sphere/plane geometry, as described in Ref. [67] (Eq. 2.4). The cone is 1 µm-high with a 10 • open angle. The sphere has a radius R (cf. Fig. 1), the value of which is specified in the following. The second vdW contribution is calculated between the surface adatom and the spherical apex of the tip 66 :

F vdW ad tip (d) = -4πC vdW ρR 3 1 (d + 2R) 4 d 3 + 1 (d + 2R) 3 d 4 , (4)
with d being the distance between the sphere' surface and the adatom center measured along the axis connecting it to the sphere's center. When assuming an interface consisting of similar materials, for example a gold adatom interacting with a gold tip, C vdW scales as α 2 and is connected to the Hamaker constant of the tip-surface interface, H, according to

H = π 2 ρ 2 C vdW ,
where ρ stands for the atomic density of the material. For gold, ρ Au = 6 × 10 28 at.m -3 and for metals, it is known 66 that H = 1 -5 × 10 -19 J. In the following, we take H Au = 3 × 10 -19 J. Therefore C vdW = 0.8 × 10 -77 J.m 6 .

The capacitive force involving the capacitance of the void capacitor C 0 was derived from

Ref. [67] as well (Eq. ( 2.3)). The expressions of the potentials φ 0 and φ im are derived by neglecting the influence of the conical part of the tip apex. They are calculated by the process of image charges of the charge q between the counter-electrode and the sphere standing for the apex of the tip, similarly to the method extensively detailed in Refs. [70-72].

D. Short-range part of the total force

The SR part of the total force is:

F SR = F vdW SR + F el SR = F vdW SR - ∂ ∂z - 1 2 -→ p tip • -→ E ext tip - 1 2 -→ p ad • -→ E ext ad (5)
F vdW SR depicts the SR London dispersion force between the interfacial atoms, that is the tip atom and the adatom. The latter component occurs whether the adatom is charged or not.

Most of the time, when the interacting atoms may be assumed as punctual objects, the interaction is described by a Lennard-Jones potential scaling as r -6 . Hence: F vdW SR = -6C vdW /r 7 . According to the estimated value for C vdW , one has used: F vdW SR (nN ) -500/r 7 ( Å). We now focus on the description of the dipolar force. The expressions of the microscopic dipoles -→ p tip and -→ p ad and of the external fields -→ E ext tip and

-→ E ext ad are to be found. However, the microscopic dipoles arise from the local electric field felt by the tip atom or by the adatom [START_REF] Kittel | Introduction to Solid State Physics[END_REF] including, in addition to the external field, the mutual influence of the dipoles on each other, thus ensuring the self-consistency of the model. Therefore,

- → p tip,ad = 0 α -→ E loc tip,ad , (6)
with

-→ E loc tip,ad = -→ E ext tip,ad + -→ E dip tip,ad . (7)
-→

E dip tip ( -→ E dip ad )
is the dipolar field induced by the polarized adatom (the polarized tip atom), its set of image dipoles, but also by the image dipoles of the polarized tip atom (the polarized adatom), influencing the tip (the adatom). The preceding equation states that the polarizabilities of the tip atom and of the adatom are similar, consistent with our assumption.

Physically, the external fields

-→ E ext tip,ad are both proportional to V b and q [visible in Eqs. (A.13), (A.14), (A.15) and (A. 19), (A.20) in the Appendix). Besides, the self-consistency of the approach implies that -→ p tip,ad stem from polarization effects, that is from a field, or a potential, able to trigger them. In turn, the dipolar fields

-→ E dip tip,ad = -→ 0 if and only if -→ p tip,ad = -→ 0 .
This requires either V b = 0, or q = 0 (in this case, the Coulombic potential created by the charge is = 0). For instance, the case V b = 0 and q = 0 gives rise to dipolar forces because a neutral adatom may be polarized. Conversely, the peculiar case V b = 0 and The microscopic origin of the dipolar forces described here differs significantly from the origin of vdW forces between atoms or molecules, among which are London dispersion forces (generic term F vdW in the model), but also Debye or Keesom forces (not described as they are irrelevant with the present setup). London dispersion forces stem from the net dipolar interaction resulting from the instantaneous electronic polarization of the interacting atoms.

q = 0, implies -→ E ext tip,
They occur whether the interacting atoms are charged or not and whether a bias voltage is applied or not. Debye forces stem from the interaction between an atom or a molecule carrying a permanent dipole and an induced dipole. Keesom forces stem from the interaction between two permanent dipoles subject to random orientation due to thermal motion. Here, dipolar forces do rely on the electric biasing of the electrodes and on the presence of net charges within the tip-surface interface because the induced fields polarize the interfacial atoms that subsequently self-influence themselves. Furthermore, they do not require the atoms to carry permanent dipoles and are not sensitive to thermal effects because the atomic positions are considered to be frozen and the electronic polarization is independent of the temperature.

III. ANALYTICAL RESULTS AND COMPARISON WITH EXPERIMENTAL

RESULTS

A. SR-electrostatic versus SR-vdW forces: strength, lateral resolution and distance dependence

In Fig. 2(a), the SR force including vdW and dipolar contributions, F SR = F vdW SR + F el SR , is drawn as a function of the lateral position of the tip depicted by its angle θ with respect to the charged adatom (θ = 0, cf. Fig. 1). The force curves have been computed with V b = 0 for two tip-surface separations, z = 4.8 and 5.1 Å and above q = -1 (charged adatom) and q = 0 (neutral adatom). The unchanged parameters are: polarizability for gold α = 6.78 Å3 , dielectric thickness h d = 3.0 Å (two NaCl layers), and ef f = 1.4 0 . These force profiles are to be compared to the experimental ones shown in Fig. 2(B) of Ref. [56] measured on Au - and Au 0 , for which the LR background force has been subtracted. In addition, we have assumed the two atoms are far enough apart ( 27 Å in the experimental work) to avoid a potential cross-influence in the resulting force measured on top of each of them.

Let us consider the case q = 0 first [neutral adatom, green curves in Fig. 2 , ont top of q = -1 (black curves) and q = 0 (green curves). The full width at half-maximum for the black curve ( q = -1, z = 4.8 Å) is ∆θ = ±0.44 rad. The comparison between the numerical results and the experimental curves from Gross et al. 56 is good (cf. text). (b) Distance dependence of the tip-adatom SR force computed with V b = 0 on top of q = -1 (black dots) and on top of q = 0 (squares). The latter curve merely consists of the SR-vdW force component (∝ z -7 ). The fit of the curve measured on top of q = -1 (continuous red curve) has been obtained with a z estimate of the force gives -8.5 pN whereas the reported experimental force rather is -14 pN. This discrepancy may be reduced if one considers a sphere-sphere vdW interaction instead of vdW interaction between atoms treated as punctual objects. The vdW interaction force between spheres may be derived from Eq. (1) in Ref. [76]. With two gold atoms (144 pm radius) 4.8 Å far apart then F vdW,max SR =-12 pN. Even though the calculations look more accurate, we keep on using the vdW interaction between punctual objects in the following in order to: (i) be consistent with the punctual-like description of the charged adatom for the electrostatic interaction and (ii) not bring further uncertainty in the calculations, as the value of the Au -radius was not found.

When now considering the curve computed on top of q = -1 at z = 4.8 Å [continuous black curve in Fig. 2(a)], the numerical estimate gives -24.5 pN and the experimental value of the force reaches -27 pN. At z = 5.1 Å(dotted black curve), the force has decayed down to -16.9 pN, corresponding to a max/min ratio of 1.45 within ∆z = 30 pm only. Therefore in the present situation, SR-vdW forces are not able to account for the experimental variations observed on top of Au -, which states that SRE forces weight significantly in the total force upon the polarizability and the charge state of the interfacial atoms.

In Fig. 2(a), it is restated that the calculations are performed with (i) both q = 0 or q = -1 adatoms lying on the dielectric' surface and hence located at the same separation with regard to the tip atom when placed above them and (ii) with the same polarizability α for them. Thus, this allows us to strictly compare SRE and SR-vdW forces, while getting rid of the topographical influence and of the chemical nature of the underlying adatoms, respectively. We are aware that these assumptions are problematic. Indeed, regarding issue (i), it is known that, when adsorbed on a NaCl thin film, Au -embeds within the surface by about 0.4 Å, whereas the neutral gold atom remains at the surface 77 . Regarding issue (ii), it is obviously expected α(Au -) > α(Au 0 ). Nevertheless, a rather good agreement is observed between the calculations and the experiments. We think this stems from the two following effects that actually balance each other: the underestimation of the tip-Au -separation, which implies an overestimation of the force [as it decays quickly with z; cf., Figs. 2(a of the dc bias. For that purpose, the expression of ∆f is derived from the force F according to a simplified expression, only valid when using small oscillation amplitudes 78,79 as this is Note also that the way the LCPD is shifted (larger values with negatively charged adatoms and vice-versa) is consistent with the experimental findings. To assess the influence of SRE forces, the curves have been plotted with (continuous lines) and without (dotted lines) the SR components.

SRE forces do not contribute significantly to shift the LCPD upon the charge state of the adatom (similar positions of the parabolas' maxima), but rather shift the total force. Therefore, with the present setup, the LR-electrostatic interaction is solely responsible for the observed LCPD shift (∆V = +80 or -100 mV). (c) ∆f (V b ) curves computed above q = -1 (black curve) and q = +1 (red curve) upon similar conditions as (b), except that the LR components are set equal to 0.

The shaded area depicts the range of bias shown in (b). Over that range, the variation of the SRE-connected ∆f is very small, despite a strong splitting between parabolas' maxima.

the case in Gross' experiments:

∆f = - f 0 2k ∂F ∂z . (8)
In Fig. 3, ∆f vs. V DC curves are reported. They have been computed on top of an adatom with q = +1, q = 0 or q = -1. We have used the experimental parameters given in

Ref. [56]: k = 1800 N/m and f 0 = 23000 Hz. The curves are computed at constant height, z = 5.8 Å with a tip radius R = 50 Å and an additional adatom with a charge state q = -1, located 10 Å aside the one on top of which the spectroscopic measurement is performed [cf. Fig. 3(a)]. This is consistent with the geometry of the spectrometric measurements performed by Gross shown in Figs. 3(A),(B) and (C). Note that the spectroscopic curve computed on top of q = +1 is for predictive purpose as no equivalent curve is reported in the experimental work. Let us focus on the curves corresponding to the total force first [continuous lines labeled as LR+SR vdW +SRE in Fig. 3(b)]. The qualitative behavior between numerical and experimental results is good and similar trends are observed. First, over a similar range of voltage (±400 mV), the numerical values of ∆f range from -2.4 down to -4.4 Hz, while the experimental range is -2.4..-3.3 Hz. We assign this discrepancy to our assumption of a sphere-plane capacitor for the calculation of the potential φ 0 and φ im in the LR-electrostatic interaction. Second, negatively charged adatoms have a larger LCPD than neutral and positively (not shown in the experimental results) charged ones.

The gap voltage between q = 0 and q = -1 is ∆V =+80 mV, to be compared to the experimental value: ∆V +30 mV. Third, a larger attractive force occurs on top of q = ±1 (more negative shift) than on q = 0, consistent with the experimental findings and with the results of the preceding section. Fourth, the maxima of the numerical curve measured on top of q = 0 is not centered around 0, but rather around V DC = +80 mV, to be compared to the 180 mV-experimental value. This effect stems from the electrostatic influence of the additional adatom q = -1 located nearby in the capacitor, which contributes to the positive shift of the LCPD, consistent with what is observed on top of q = -1. This is also the effect that makes the LCPD calculated on top of q = +1 almost zero since the total charge within the capacitor is then zero.

We now discuss the relative weight between LR and SR forces in ∆f (V b ) spectroscopic curves. To first order, dipolar forces scale as V 2 b (cf. Appendix), but the bias-dependence of the total force as well stems from the LR-electrostatic part of the interaction. In order to understand what force component is the most significant in the bias dependence, spectroscopic curves computed without SR forces, that is only with the LR components, are drawn in the figure [dotted curves labeled as LR in Fig. 3(b)]. The parabolas are not bias-shifted, but ∆f -shifted (less negative values), upon removal of SRE forces. Indeed, when investigating the bias dependence of the mere SRE forces above q = -1 and q = +1 upon similar conditions [geometry, height; cf. Fig. 3(c)], it is found that over the range ±400 mV (shaded area in the figure), the curvature of the parabolas is very weak, corresponding to ∆f changes of about 0.1 Hz. However, a large splitting between parabolas' maxima (∆V DC = 19 V) is observed, which states that SRE forces are extremely sensitive to the charge state of the adatom as soon as LR-electrostatic forces are not present in the tip/dielectric/counterelectrode capacitor. Conversely, when these are present, the spectroscopic curves become weakly influenced by SRE components, which fully agrees with the pioneer work by Terris et al. 57 . In this case, any attempt to accurately fit ∆f (V b ) experimental data requires a precise calculation of the electric field within the capacitor 40,80,81 , and therefore a precise knowledge of its geometry.

C. SR-versus LR-electrostatic forces: influence on the lateral resolution

We here address the issue of the lateral resolution as a function of SR-and LR-contributions. The conditions to achieve high-resolution imaging rely on SR interactions, including dipolar as well as SR-vdW forces, but also on the tip sharpness, that is the size of the tip apex.

For that purpose, we mimic the experimental situation shown in Fig. 2(A) of the work by Gross et al. A negatively charged adatom ( q = -1) lying on the thin dielectric is placed at a distance d = 25 Å from a neutral adatom ( q = 0). In Fig. 4, ∆f profiles are reported as a function of the tip size and force components. These are derived from the corresponding force profiles calculated at z = 4.8 Å with V b = 0 according to the same method as before.

For each case reported in the figure, we have plotted the total force including LR and SR components (red curves) and solely the LR component (including vdW and electrostatics, black curves). The LR force is calculated for tip apexes corresponding to: (i) a hemispherical apex with R=20 Å [Fig. 4 Let us first discuss the influence of the LR interaction on the calculated ∆f profiles. As expected, the LR interaction increases with the size of the macroscopic tip and produces a ∆f offset ranging from -1 Hz in Fig. 4(a) up to -47 Hz in Fig. 4(d). At the distances used in the experimental setup, it is the LR electrostatic interaction which dominates over the LR-vdW one for the sharpest tip. This can be seen in Fig. 4(a) where the different contrast of the black curves above Au 0 and Au -is caused by the LR electrostatic interaction only.

In the present case of zero bias (V b = 0, and hence φ 0 = 0) and a tip radius R which is much larger than the distance z (R > 5z), the LR electrostatic interaction only depends on the charge and the induced potentials φ im [see Eq. ( 2)] and preserves a more or less constant strength ( 2 Hz) for all tip sizes. It is always centered around the charged adatom, however, its lateral confinement is decreasing with increasing tip size. Other than the LR electrostatic forces which are constant in the present setup, the LR-vdW forces scale with the volume of the tip. This is why for the flat tip calculated in Fig. 4(d) it is the LR-vdW which produces the largest part of the interaction. The LR-vdW is not sensitive to the charge; consequently, the maximum contrast is now centered between the two adatoms and no longer above the charged one. Please note that the LR-vdW part contains two contributions. One describes the interaction between the adatoms and the macroscopic tip; the other describes the interaction between the surface and the macroscopic tip. It is the first one which causes the parabolic shape of the black curve in Fig. 4(d), the latter one just adds a constant offset to the calculated ∆f .

All the previously described LR interactions produce a more or less structured background in the ∆f profiles due to a convolution with the macroscopic tip shape. As soon as SRE forces are present, the convolution with the tip apex does not prohibit resolving charged adatoms. The main reason is that SRE forces are able to spatially and selectively pin the force on top of each adatom because they prevail on SR-vdW forces, as discussed in subsection II D. However, the lateral resolution depends on the exact shape of the apex.

Surprisingly, a tip with an apex the size of which is significantly larger than the separation between the species to be resolved, provides the maximum of lateral resolution. Indeed, the influence of LR components is then averaged over the area of the tip, whereas the spatial distribution of the SR dipolar and vdW forces remains unchanged as they are not apexsize dependent. Conversely, this tip will provide no sensitivity to the charge state of the adatoms because the LR-electrostatic contribution will weight too much in the resulting bias dependence of the total force, as described in subsection III B. A very sharp tip [case (i)] will provide high resolution as well, as the one used by Gross et al. to perform their experiments, R 5 -10 Å. Between these extreme cases, the resolution will depend on the size of the tip apex, as discussed, for instance, in Ref. [82].

IV. SUMMARY AND CONCLUSION

We present a model which explains in detail the nc-AFM experiments reported by Gross et al. 56 . In these experiments, nc-AFM is used first to image the topography of charged and neutral gold atoms adsorbed on a thin dielectric layer deposited on a metal substrate.

Second, frequency shift ∆f as a function of bias voltage V b curves are used to identify the charge state of these adatoms.

Our model mimics the experimental tip-sample geometry and accounts for both LR-and substrate covered by two monolayers (3 Å) of a dielectric film. On both, tip and surface, there These findings allow us to draw conclusions regarding experimental strategies to improve the spatial resolution on isolated charged species and/or high spectroscopic sensitivity to their charge state. In particular, the proposed experimental geometries might be the key to success in order to investigate surfaces that will be influenced or even destroyed when the tip apex approaches within the range of atomic binding distances.

• The topographical appearance of charged adatoms adsorbed on a dielectric relies on SRE forces, their high-resolution imaging requires a polarizable tip apex, but not mandatorily ultra-sharp. The latter issue however, would obviously facilitate the high-resolution. Polarizable tips might be achieved by using a properly-grafted polarizable molecule at the tip apex, or by using a tip coating consisting of a high-K material, for instance.

• The charge sensitivity in either the spectroscopic signature or KPFM measurements is hindered by LR-capacitive forces. In order to increase sensitivity, SRE contribution should be predominant. Thereto, two suggestions can be made. The first one is to increase the SRE contribution by using both, a polarizable tip apex as well as a sharp tip in order to increase the local electric field. The second concept is to make the tip/counter-electrode capacitance as small as possible. For that purpose, using a bulk insulator instead of a thin layer might be advantageous.

Choosing such experimental conditions would allow, for example, the detection of the local charge distribution within a single dipolar molecule.

FIG. 5: Geometry of the parameters defined in the model and description of the image dipole process for the calculation of the dipolar fields, E dip ad and E dip tip . For that calculation, the metallic tip apex is assumed as being planar (cf. text). Owing to the image process, this plane as well as the counter-electrode plane are antisymmetrical. Indexes (0) refer to the real dipoles. The set of image dipoles within the tip apex are numbered with a prime ( ). For the calculation of the dipolar fields, we have solely used the set of images dipoles depicted in the figure (cf. text).

where

-→ E ext tip and -→ E ext ad are the external electric fields influencing the foremost tip atom and the surface adatom, respectively. These are produced by the biased void capacitor, the charged adatom and its set of image charges within the metallic electrodes.

In addition to the latter external fields, the self-consistency of the problem is ensured by the assumption that the adatom is as well polarized by the dipolar field of the tip atom and vice versa:

- → p tip,ad = 0 α -→ E loc tip,ad , (A.10) with:
-

→ E loc tip,ad = -→ E ext tip,ad + -→ E dip tip,ad (A.11)
-→ E dip tip is the dipolar field polarizing the foremost tip atom. It originates from the following:

• the polarized adatom

• the set of image dipoles of the polarized adatom within the metallic electrodes

• the set of image dipoles of the polarized tip atom within the metallic electrodes -→ E dip ad is the dipolar field polarizing the surface adatom. It originates from:

• the polarized tip atom

• the set of image dipoles of the polarized adatom within the metallic electrodes

• the set of image dipoles of the polarized tip atom within the metallic electrodes

The dipolar field in a point M of the space produced by an isolated dipole -→ p located in a point O, such that --→ OM = r e r , is:

-

→ E dip (M ) = 1 4π 0 3( -→ p • e r ) e r --→ p r 3 (A.12)

Once

-→ E ext tip,ad are calculated, the combination of Eqs. (A.12), (A.11) and (A.10), allows the derivation of a linear set of equations for the x and z components of -→ p tip and -→ p ad .

Because SRE forces have a very narrow lateral extension, the calculation of the dipolar fields is performed by assuming a planar geometry of the tip-surface interface, which means we assume R z. Hence, due to the effect of the image dipoles process, the planar electrodes become antisymmetrical planes. Since -→ p tip belongs to the tip apex plane, it is oriented perpendicularly to this plane, that is, along the z axis. As discussed in the text, we group -→ p tip and its image with regard to the tip plane, which ultimately implies to replace p tip,z with 2p tip,z in the calculations. Conversely, -→ p ad cannot be projected along a single direction.

Calculating

-→ E dip tip and

-→ E dip ad by using the image dipoles process leads formally to an infinity of image dipoles of -→ p tip and -→ p ad . Nevertheless, because the dipolar field decreases as r -3 [Eq. (A.12)] and because the distance between each image dipole and the real dipole increases as the number of images under consideration increases, we restrict the calculation of the dipolar field to the image dipoles depicted in Fig. 5. This configuration stands for the minimum set of images that ensures the antisymmetry condition on the tip apex as well as the coupling between -→ p tip and -→ p ad and their first images. Thus, with this level of approximation,

-→ E dip tip is derived by only considering interactions with -→ p When assuming that the void capacitor has a planar geometry, the expression of the electric field becomes trivial:

-→ E ext,0 = V b h d + z z. (A.14)
The electric field produced by the image charges of the charge may be derived by recurrence,

-→ E ext im ad = q 4π ef f ∞ n=0 1 [2z + 2d n ] 2 - 1 [2h d + 2d n ] 2 z, (A.15)
with d n = n(h d + z) and n ∈ N. Therefore, -→ E ext ad has a single vertical component. When considering both the convergence of the series and the computing time for the calculations, it is found that n = 20 is an optimal choice for shrinking the series.

(ii)

-→ E dip ad states for the dipolar field influencing the surface adatom. The vectorial projection along x and z axes gives: The tip dipole has solely a vertical component, as expected from the symmetry conditions.

The set of linear equations (A.17

  [47], Sadewasser et al. have shown that polarization effects were responsible for the atomic scale KPFM contrast on the Si(111)-7×7 including Pb adatoms, as well. Their results point toward a relation between the changes in the LCPD over distinct atomic sites and the variations of the surface local electronic structure due to a charge polarization induced by the tip-surface interaction. More recently, Bieletzki et al.54 have also stated that the strong dipole in MgO thin films deposited on Ag(111) leads to specific KPFM and topographical imaging conditions depending on the charge state of the tip. More recently even, Masago et al.55 , using a theoretical approach based on DFT calculations, stressed the influence of the polarizability of atomic species in the apparent LCPD on Si surfaces and reach markedly similar conclusions as those of the present work.Most of the aforementioned results indicate that polarization effects play a key role in the occurrence of SRE forces with two consequences in nc-AFM and KPFM: (i) They contribute to the total interaction force and are thus involved in the ∆f signal used to perform "topographical images" and (ii) despite about a few tens of pico-newtons43,56 , they must be considered for interpreting the LCPD atomic scale contrast owing to their cross-coupled distance-and bias voltage-dependence. To a first approximation, SRE forces should be proportional to the square of the applied bias voltage V b . Thus, in both KPFM and spectroscopy curves where the frequency detuning ∆f is recorded as a function of the applied bias voltage V b [∆f (V b ) parabolas], it should be possible to detect the SRE contributions to the total interaction. Besides imaging by means of KPFM, the aforementioned single point spectroscopy ∆f (V b ) parabolas and similar detection methods are more frequently used to detect local charges and charge transfer. Since the first experiments by Terris et al.

  -→ p tip,ad = -→ 0 . Then no dipolar forces occur, and the resulting SR force is only driven by the SR-vdW interaction.The formal expressions of the external and local electric fields are given in appendix. The lateral (x axis) and vertical (z axis) expressions of -→ p ad and -→ p tip are derived by solving a set of linear equations [cf. Appendix, Eqs. (A.17) and (A.22)]. Finally, the vertical component of the SRE force, F el SR , may be derived by differentiation with respect to z [Eq. (5)], and thus the expression of the total force may be obtained. The problem thus rationalized relies on two main parameters, namely, the polarizability of the interfacial atoms α and the Hamaker constant of the interface H.

  (a)]. No dipolar force occurs because V b = 0 and q = 0. The resulting interaction is only driven by the SR-vdW force, as stated in Sec. II D. We stick to the experimental values of the force measured at z = 4.8 Å in Ref.[56] [yellow curve in Fig. 2(B) of the reference]. The numerical

FIG. 2 :

 2 FIG. 2: (Color online). (a) Tip-adatom SR force including SR-vdW and dipolar contributions as a function of the angular position of the tip with respect to the adatom located in θ = 0. The force profiles have been computed with V b = 0 and z = 4.8 Å (continuous curves) and 5.1 Å (dotted curves), ont top of q = -1 (black curves) and q = 0 (green curves). The full width at half-maximum

- 6 function.

 6 The experimental values of the force derived from Gross et al. are depicted with blue squares. The horizontal error bars depict a possible uncertainty in the z position of the tip of 2% only. The discrepancy between SRE and SR-vdW forces is well visible.

) and 2 (

 2 b)], and an underestimation of the force on top of Au -, as all the SR forces increase with α.The force profile as a function of the lateral position of the adatom allows us to draw conclusions regarding the lateral resolution of neutral and charged adatoms in the experimental images. This requires first to define a criterion stating unambiguously what the lateral resolution is. We stick to the following criterion: two objects separated by a distance d over the surface are assumed as being resolved by the tip if the full width at half maximum (FWHM) of their respective force profile [∆x , respectively] do not overlap, that is d > (∆x )/2. In the range 4.8-5.1 Å, the calculated FWHM force profile on the negatively charged adatom or on the neutral one does not vary much: θ F W HM ±0.44 rad. Therefore, at z = 4.8 Å, ∆x F W HM (Au -)=∆x F W HM (Au 0 )=2z × tan(θ F W HM ) = 4.5 Å. The experimental values are: ∆x F W HM (Au -) 5 Å and ∆x F W HM (Au 0 ) 7.5 Å. If the comparison to the experiments is good in the case of Au -, which testifies to the importance of the SRE forces contribution, there is a significant discrepancy with Au 0 . At this stage, this is not yet understood. But, regarding the spatial extension of the force profile on a single electron charge, the "apparent size" of a charged adatom is a few Angströms only. This is a somewhat unexpected result in that a naive view of the problem might have expected a usual Coulombic z -1 dependence, thus preventing any high-resolution imaging of these objects from being performed. In Fig. 2(b), the distance dependence of the total SR force on top of an adatom with q = -1 is reported (black dots) and compared to the experimental data from Gross et al. (blued squares). These have been derived for few tip-surface separations given in the reference [Fig. 2(B); z = 4.8, 5.15, and 6.0 Å]. The agreement between the experimental and the numerical data is still good, within however deviations due to the limitations of the model and to the experimental error bars. For instance, the z error bars shown in Fig. 2(b) correspond to ±2% only. The best fitting function for the force was found to be z -6 (continuous red line). For comparison, we have as well drawn the SR-vdW curve (i.e. on top of Au 0 ) scaling in z -7 (black squares), the behavior of which is significantly different. This trend is fully consistent with the one shown in Fig. (S1B) of the supporting online material of Ref.[56]. Thus, SRE forces decay almost as fast as SR-vdW, but their strength may prevail on the SR-vdW contribution upon the polarizability and/or the charge state of the interfacial atoms. B. SR-versus LR-electrostatic forces: influence on the charge state characterization A major aspect of the work by Gross et al. deals with the detection of single electron charges by means of the measurement of the LCPD. The authors showed that the LCPD might selectively be shifted upon the sign of the charge of the underlying adatom (positive, neutral or negative). We hereafter detail the bias voltage dependence of the total force derived from the model upon the sign of the adatom and discuss how the LR-electrostatic component may influence the detected LCPD. In order to stick to the experimental results of Gross et al. shown in Fig. (3A) of the reference, the spectroscopic curves are reported in terms of frequency shift (∆f ) as a function

FIG. 3 :

 3 FIG. 3: (Color online). (a) Geometry used to compute the spectroscopic curves shown in panel (b). (b) ∆f vs. dc bias spectroscopic curves computed for q = +1 (red), q = 0 (green), and q = -1 (black) at z = 5.8 Å. The tip apex's radius has been set equal to R = 50 Å. The numerical curves are found to be in good agreement with the experimental results measured on top of Au 0 and Au -.

  (a), R d], (ii) same as (i) with R=45 Å [Fig. 4(b), R > d], (iii) same as (i) with R=200 Å [Fig. 4(c), R d], and (iv) same as (i) with R = 1000 Å [Fig. 4(d)].

FIG. 4 :

 4 FIG. 4: (color online). ∆f profiles computed from the corresponding force profiles calculated at z = 4.8 Å with V b = 0. The charged adatom ( q = -1) is placed in x = 0 and the neutral adatom is placed at a distance d corresponding to x = -25 Å. Cases a-, b-, c-and d-correspond to a tip with hemispherical apex R=20 Å (R d), R=45 Å (R > d), R=200 Å (R d) and a flat apex with an area equivalent to a hemispherical apex of radius R = 1000 Å, respectively. For each case, the total force (red curves) or the total force minus the SR components (black curves) has been used to derive the ∆f profile.

  is a (polarizable) gold atom adsorbed. The SR chemical binding forces and displacements of atoms are not included since in the experiments, the closest tip-sample distance is larger than 4.5 Å. All material constants in our model are chosen according to values given in the literature. The model is based on the calculation of the electrostatic energy for the complete geometrical system. In addition to the well known terms of LR-and SR-vdW, as well as the LR-electrostatic terms, our model includes a less common term of SR-electrostatic (SRE) forces stemming from the self-consistent description of electronic polarization at the tipsurface interface including the contribution of image charges and image dipoles. These SRE interactions, so called dipolar forces, have a microscopic origin due to the polarizable atoms adsorbed on the tip front end and on the surface. The polarization of these atoms is due to the biased electrodes (metallic tip and counter-electrode), as well as to charges within the tip/counter-electrode capacitor.The previously described SRE interactions are weak in comparison to SR binding forces.The experiments by Gross et al. were performed at distances larger than 4.5 Å and therefore serve as model case experiments since SR binding forces are negligible. In this model case, it is found that the interaction responsible for the atomic scale topographical contrast in nc-AFM is not determined by SR-vdW, but rather dominated by SRE interactions. However, for the identification of the charge state of the adsorbed atoms by means of ∆f (V b ) spectroscopy, our calculations clearly point out that the experimentally observed variation of the LCPD is dominated by the LR-electrostatic interactions. Thus, a different interaction allows for the discrimination of charged adatoms in topography imaging and in ∆f (V b ) spectroscopy, respectively.

  electric field felt by the charged adatom: E loc ad We seek the expression of the local electric felt by the charged adatom -Eq. (A.11)].(i) The external electric field acting on the charged adatom -→ E ext ad consists of two components: (a) the field of the void capacitor that may locally be assumed as a planar capacitor (R z), referred to as -→ E ext,0 , and (b) the field produced by the image charges of the charge within the metallic electrodes of the capacitor, that is the spherical apex of the tip and the counter-electrode, referred to as -

r 3 1 ( 1 (E

 311 cos θ sin θ -6p tip,z r 3 cos θ sin θ + 2p ad,x 2x) 3 + 1 (2z) 3 , 2p tip,z r 3 (3 cos 2 θ -1) + 2p tip,z r 3 (3 cos 2 θ -1) + 2p ad,z 2h d ) 3 + 1 (2z) 3 . (A.16)Therefore, the equation giving the expression of the dipolar moment of the adatom becomes:-→ p ad = 0 α tip,z ; p ad,x ) E ext ad + E dip ad (p tip,z ; p ad,z ) (A.17)As expected, the dipole of the polarized adatom has two vectorial components. The x component may be derived and put in a simple shape, that is at lowest order in α. We get:p ad,x = α 4π 6p tip,z cos θ sin θ r 3 -cos θ sin θ r 3 (A.18)Therefore, it may be seen that the x component of the dipole of the polarized adatom stems from the z (i.e. vertical) component of the tip dipole.C. Local electric field felt by the tip atom: E loc tipWe now seek the expression of the local electric felt by the tip atom: consists of three components: (a) the field of the void capacitor -→ E ext,0 already found, (b) the field of the image charges of the charge within the metallic electrodes of the capacitor, and (c) the Coulombic field due to the charged adatom. However, the two latter components are derived in a row and give the expression of the electric field referred to astan θ) 2 + (z + 2d n ) 2 ] 3/2 -2h d + z + 2d n [(z tan θ) 2 + (2h d + z + 2d n ) 2 ] 3/2 z (A.20) with d n = n(h d + z).For the same reason as before, the calculations are performed withn = 20. Contrary to -→ E ext im ad ,-→ E ext ad+im tip is not zero for n = 0, which actually gives the Coulombic electric field produced by the charge influencing the tip.(ii)-→ E dip tip states for the dipolar field influencing the tip atom. We get: h d + z)) 3 + 6 r 3 (p ad,x cos θ sin θ + p ad,z (cos 2 θ -1)) z (A.21) Therefore, the second equation giving the expression of the dipolar moment of the tip atom is: ext,0 + E ext ad+im tip + E dip tip (p tip,z ; p ad,x ; p ad,z ) (A.22)

  ) and (A.22) gives the expressions of p ad,z and p tip,z , which ultimately gives the expression of F el SR [right-hand side in Eq. (5)].

Acknowledgments

The authors would like to thank Leo Gross and Gerhard Meyer for valuable comments and acknowledge financial support of the Centre National de la Recherche Scientifique, of the Région Provence Alpes Côtes-dAzur, as well as of the Agence Nationale de la Recherche with the PNANO project MolSiC (ANR-08-P058-36).

APPENDIX

The geometry of the parameters defined hereafter is detailed in Fig. 5.

A. Summary of the analytical method and main assumptions

The vertical component of the SR dipolar force to calculate derives from the classical equation: