Ab-initio calculation of the real contact area on the atomic scale

We present a novel approach to determine the onset of contact between a tip and a surface. The real contact area depending on the distance is calculated using Bader's quantum theory of atoms in molecules. The jump to contact, which is often observed in atomic force microscopy experiments, is used as an indicator for the initial point of contact, which in turn is defined by atomic relaxations and thus without the need of external parameters. Within our approach the contact area is estimated by evaluating the zero flux surfaces between the touching Bader-atoms, where the necessary electronic density cutoff for the Bader-partitioning is calculated to depend on the initial point of contact. Our proposed approach is therefore completely ab-initio and we are able to define and calculate the real area of contact without imposing restrictions or free parameters. As a prototype system we choose a tip made of a ten atom tungsten pyramid above a moir\'{e} layer of graphene on an fcc iridium (111) substrate. We find that the contact area depends exponentially on the effective distance between the tip apex and the surface atom directly below within the atomically relaxed nanosystem.

We present a novel approach to determine the onset of contact between a tip and a surface. The real contact area depending on the distance is calculated using Bader's quantum theory of atoms in molecules. The jump to contact, which is often observed in atomic force microscopy experiments, is used as an indicator for the initial point of contact, which in turn is defined by atomic relaxations and thus without the need of external parameters. Within our approach the contact area is estimated by evaluating the zero flux surfaces between the touching Bader-atoms, where the necessary electronic density cutoff for the Bader-partitioning is calculated to depend on the initial point of contact. Our proposed approach is therefore completely ab-initio and we are able to define and calculate the real area of contact without imposing restrictions or free parameters. As a prototype system we choose a tip made of a ten atom tungsten pyramid above a moiré layer of graphene on an fcc iridium (111) substrate. We find that the contact area depends exponentially on the effective distance between the tip apex and the surface atom directly below within the atomically relaxed nanosystem.

INTRODUCTION
Historically the introduction of the concept of a real contact area in 1939 by Bowden and Tabor, which is substantially smaller than the nominal one, was a huge step forward in understanding the laws of friction [1]. Since then a multitude of methods have been used to predict the pressure dependent real contact area for rough surfaces.
In 1881 Hertz laid the foundations of contact mechanics by describing the junctions between non-adhesive, homogeneous elastic solids of simple shapes [2]. Nearly a hundred years later Johnson, Kendall, and Roberts (JKR) included short-range adhesion forces inside the contact which lead to larger contact areas compared to Hertz's model [3]. In contrast, Derjaguin, Muller and Toporov (DMT) assumed that the Hertzian contact area remains undeformed while long-range adhesion forces are acting outside the contact zone [4]. Tabor was able to show that both the JKR and the DMT model can be viewed as limiting cases of a more general model, with JKR suitable for soft materials with large adhesion and DMT for hard materials with low adhesion [5]. Additional work on this unification was done by Maugis using Dugdale potentials [6]. The analytical solution by Maugis, however, produces rather cumbersome equations, which can be approximated with high accuracy by the generalized transition equation derived by Carpick, Ogletree and Salmeron [7]. Their result also describes the transition between the DMT and the JKR models, but with simpler expressions which can be applied more straightforwardly to experimental data and only differs from the Maugis-Dugdale model within the unstable low load region. Generally these theories agree that the contact area A between a single sphere and a flat surface is a sublinear function of the load L, see e.g. the original Hertz prediction of A(L) ∼ L 2 / 3 . These continuum mechanics models are undoubtedly very successful and play an important role in both theoretical and experimental work in tribology. However, on the atomic scale, as tested for example in AFM experiments, the size of the the contact approaches the size of the involved atoms and thus models based on continuum mechanics are hardly applicable [8]. The apparent success of the Maugis-Dugdale model, which is still widely used to interpret atomic scale AFM experiments, can be attributed to its flexibility provided by three fitting parameters [9]. Hence, atomistic methods are of great value in assessing the validity of the results and provide a much needed tool to increase the understanding of the real contact area.
In recent studies and discussions a need for characterizing the contact area on an atom-by-atom basis was expressed and various strategies to estimate the number of atoms in contact were proposed [8,[10][11][12][13]. However, it is not trivial to decide at which distance two atoms are in contact or how big the resulting contact area should be. One possibility is to define a certain inter-atomic distance a 0 below which contact should be established, however, this just shifts the problem to find the correct distance a 0 . A common method is to identify a 0 , and thus the onset of contact, as the beginning of repulsion between the two observed atoms [11,14]. To this end, classical molecular dynamics (MD) studies using Lennard-Jones potentials are often employed, sometimes without the attractive part if the surfaces of interest are non-adhesive. In our view, this method is not ideal for describing the onset of contact on the nanoscale for the following two reasons: (i) if only repulsive interaction is a sign of contact, the atoms in a solid or molecule at 0 K are not in contact with each other; (ii) in AFM experiments one often observes a "jump to contact", where either the tip, or some part of the surface below or both jump towards arXiv:1505.07246v1 [cond-mat.mtrl-sci] 27 May 2015 each other because of strong attractive interactions. If now the tip support is lowered further, the distance between the surface and the tip apex might get even smaller as the chemical binding becomes stronger. It seems unreasonable to argue that the tip is not in contact with the surface after the jump, just because the interaction is still attractive.
In a study from 2009, Mo, Turner, and Szlufarska examined the contact area between hydrogenated amorphous carbon tips of up to 30 nm radius and a flat hydrogenated diamond surface [10]. For this large scale finite temperature (300 K) MD study they used a reactive empirical bond-order (REBO) potential [15] to model the chemical forces and added an analytical switching function to include van der Waals (vdW) like long-range forces. The multi-asperity picture of nanoscale contact presented in their publication relies on the assumption that contact is established between the atoms that are interacting chemically through the REBO potential while a much larger part of the tip is attracted to the surface via the vdW forces. These ideas are further discussed in a follow-up publication by Mo and Szlufarska [16]. Now an atomic contact area A at is attributed to every chemically interacting atom of the tip leading to the total real contact area A real = N at A at , with N at being the number of involved atoms. Due to atomic scale roughness in the amorphous tip this real contact area may be significantly smaller than the expected contact area of a smooth asperity, A asp , which is defined as the envelope over the contact points. Concluding, Mo, Turner and, Szlufarska pointed out how important atomic corrugations are for an accurate estimation of the real contact area, and they underlined the need for accurate computational approaches. Nevertheless, a few points remain to be analyzed further. In a realistic, continuous potential it might be hard to distinguish between long-range dispersion forces and chemical forces, thus making it difficult to define the number of atoms in contact. Furthermore, the inherent assumption that all contributions from single atomic contacts A at are of equal size might not hold in systems with varying local load.
In the following, we will introduce an ab-initio approach to estimate the distance-dependent real contact area between a tip and a surface. To prove the feasibility of our approach, we choose a realistic and thus rather complex system, which has previously been used to explain contrast inversion between constant current scanning tunneling microscope and constant frequency AFM images of graphene moiré on metals [17,18]. We employ van der Waals corrected density functional calculations and use Bader-partitioning of the electronic charge density to identify accurate atomic volumes and surfaces in different chemical surroundings.

COMPUTATIONAL METHODS
The simulation cell consists of 4 layers of a 9 × 9 iridium(111) substrate covered by a 10 × 10 graphene layer forming a moiré structure with an average separation of 3.42Å and a corrugation of 0.35Å. The lattice constants obtained with the optB86b-vdW functional [19,20], a Ir = 2.735Å and a Gr = 2.465Å, are close to the experimental values of 2.71Å and 2.46Å, respectively. The mismatch of the structure is very small at 10a Gr − 9a Ir = 0.015Å [18]. The tungsten AFM tip was modeled as a ten atom pyramid with one atom at the apex, four in the next layer and five at the top. The contact site studied was an on-top position in a top-hcp region of the moiré structure. This means that the tip apex atom is positioned directly vertical over a carbon atom in a region where each carbon atom is either directly over an iridium atom or over an iridium hcp position. The simulation cell, containing 534 atoms, is shown in figure 1. Relaxations were allowed for the graphene layer and the bottom five atoms of the tip, keeping the iridium substrate and the top layer of the tip rigid at their initial relaxed positions. Relaxations of the iridium substrate during movement of the tip have been neglected due to the small binding energy (∼80 meV per carbon atom) of the mainly physisorbed graphene layer, which makes any effect of the iridium substrate on the relaxations of the graphene unlikely. The topmost layer of the tungsten pyramid, on the other hand, needs to be held fixed to control the distance between the tip and the graphene sheet.
All calculations were performed within density functional theory (DFT) employing the Vienna Ab-Initio Simulation Package VASP [21][22][23][24] using the Projector Augmented-Wave (PAW) method [25,26]. To include van der Waals (vdW) forces, which are relevant in this system, the optB86b-vdW functional was employed [19,20]. This vdW density functional has been applied to a wide range of materials and proven to be of good accuracy [27][28][29][30][31][32][33]. The Brillouin zone sampling was performed on a Γ-centered 3 × 3 × 1 k-grid, with a smearing of 0.1 eV using the method of Methfessel and Paxton to first order [34]. To ensure good accuracy for the Bader-partitioning scheme, the electronic charge density was calculated on a dense mesh of 432 × 432 × 448 points in the simulation cell. Electronic energies were converged to 10 −6 eV and the ionic relaxations were stopped after converging forces between the interacting atoms to better than 0.01 eV/Å. Following the investigation by Garhofer, who also provided us with initial structural data [18], we choose a plane wave cutoff of 300 eV, which is the minimum recommended value for carbon, and at the same time larger than the suggested value for iridium and tungsten.
To partition the electronic charge density ρ in our sim- ulation cell into single atoms we use Bader's quantum theory of atoms in molecules (QTAIM) [35]. In contrast to other similar approaches [36][37][38][39][40][41], Bader's method produces non-overlapping atomic domains with well-defined boundaries, which are perfectly suited to analyze the contact between two adjacent bodies. The necessary and sufficient condition that needs to be fulfilled to define the boundaries of a selected atom according to the QTAIM is formulated using the basic quantity in DFT, namely the electronic charge density ρ and is given as [42], ∇ρ(r) · n(r) = 0 ∀r ∈ S(r) .
Here S is the boundary surface of the atom and n(r) is the unit vector normal to this surface. The condition states that the flux of the gradient field of the charge density, ∇ρ, through the boundary surface S must vanish, which is thus called a zero flux surface. The Bader analysis in this work was performed with the code developed by Henkelman, Sanville, and Tang which is directly compatible with the format of VASP -output files [43][44][45].

RESULTS
If one seeks to define the real contact area on an atomic scale, it is quite natural to think about the size, shape, and deformations of the involved atoms. Once the size and shape of all atoms in the contact region are determined, the calculation of the real contact area is reduced to a simple summation of the regions that are in contact, provided one can distinguish unambiguously between the two contacting bodies. Since ρ formally is non-zero everywhere, the Bader-atoms at the surfaces of the contacting bodies extend into the vacuum region to infinity or until they encounter another atom. This would mean that contact between two bodies is established at all distances, which is a clearly unphysical result, unless one defines a density cutoff. This density cutoff ρ cut cannot be chosen arbitrarily, since it directly influences the contact area.
A possibility to extract a value for ρ cut is to analyze the interaction potential between the tip and the surface, divide it into a long-and a short-range part and define the contact at the onset of the short-range interaction, analogous to Mo, Turner, and Szlufarska [10]. This procedure defines a density cutoff ρ cut so that the Bader-partitioning yields contact only after the onset of short-range interactions. However, as long-range interactions are included implicitly in the exchange-correlation potential that we use (see section Computational Methods), the separation into a long-and a short-range part is not straightforward. A more unambiguous way to define the onset of contact is to analyze the atomic relaxations that happen if the tip is lowered towards the surface. We distinguish the distance for the static, unrelaxed system (d s ) and the relaxed distance (d r ), which are both measured between the tip apex atom and the carbon atom directly beneath it. For large distances no relaxations will happen, although there might be attraction due to vdW forces, and the distance d r in the relaxed system will be equal to the (static) distance d s measured before relaxing the system. At some point during the approach of the tip stronger forces will cause relaxations, which will result either in a "snap" or "jump" to contact (a phenomenon often observed in AFM experiments) if the interaction is attractive, or in a depression of the surface layer if the interaction is purely repulsive. A sketch of this process is given in figure 2. In any case, the d r versus d s curve will have a discontinuity at some distinct distance where the system begins to strongly interact. Below this distance the system will try to hold the ideal distance between tip and surface. It is straightforward to identify this discontinuity as the onset of contact.
In the examined system the interaction between the tip and the surface is attractive at the onset of contact and a jump to contact occurs between d s = 3.65Å and d s = 3.53Å (see figure 3a), such that d r is changing from 3.5Å to 2.7Å, accordingly. Most of the movement is done by the surface, which jumps upwards to meet the tip (see figure 3b). This means that when the tip support is lowered by only 0.12Å, the distance between the tip apex and the surface is reduced by 0.8Å. This allows us to define the onset of contact at d s 3.6Å and to tune the value of the density cutoff ρ cut accordingly. Note that an infinitely stiff cantilever is assumed in our calculations, as the uppermost atoms of our tip are kept atomic relaxations will not have an effect on the distance between the tip and the surface. If the tip gets closer, the surface will interact with the tip and will either jump towards the tip, if the interaction is attractive (left panel), or will get depressed, if the interaction is repulsive (right panel). This effect can be used to define the onset of contact.
rigid. In principle, the influence of a more compliant AFM apparatus on the initial jump to contact, however, could also be modeled by using a multiscale approach.
Once ρ cut is selected the determination of the real contact area for each distance is straightforward. Since the partitioning code produces only Bader volumes rather than zero flux surfaces, we have to construct the contact area from these data. To this end the Bader-volumes of both contacting bodies are added up and by pairwise comparison of the respective values for neighboring grid points a point cloud forming the contact area is generated. The contact area can now be obtained by triangulation.
We calculated contact areas for cutoff densities from 10 −3 e/Å 3 , which is the default cutoff in the partitioning code by Henkelman, Sanville and Tang [43][44][45], up to a cutoff of 10 −1 e/Å 3 . While the default value of ρ cut = 10 −3 e/Å 3 and other low cutoff densities are giving sizable contact areas for all distances, we approach the desired effect of establishing contact only for d s < 3.6Å, for a value of ρ cut ∼ 10 −2 e/Å 3 . Obviously, a very high ρ cut is unphysical, as the number of electrons that are "lost" into the vacuum region increases with rising cutoff. This means that we want to select a value that is high enough to guarantee that the contact area A is only non-zero after the snap to contact has occurred, but is otherwise as low as possible. We analyzed several values of ρ cut ranging from 7.5×10 −2 e/Å 3 to 1.0×10 −2 e/Å 3 in order to find the lowest value that still satisfies these conditions, resulting in an optimal value of ρ cut = 5 × 10 −2 electrons perÅ 3 .
Of course, changing the contact site of the tip away from a position directly above a carbon atom or to another section of the moiré structure could conceivably change the exact point of the jump to contact and thus modify the value of ρ cut . However, the obtained cutoff density of ρ cut = 5 × 10 −2 is low enough to result in an essentially flat graphene surface and hence changing the tip position should only marginally alter the computed contact area. As this paper is mainly concerned with the presentation of a new approach in defining and calculating the real contact area ab-initio, and given the rather large computational effort [46], we decided against repeating our calculations on different contact sites and calculating an average.
In a preliminary calculation of the same tip on an fcc copper (111) surface we found an optimal charge density cutoff value of 5.3 × 10 −2 electrons perÅ 3 , which is approximately the same as for the graphene/Ir(111) system.
As the optimized cutoff of 5 × 10 −2 e/Å 3 is 50 times larger than the default value it is important to check if it is still a reasonable number and does not give any unphysical results. We therefore calculated the nominal number of electrons that are assigned to the vacuum region and thus are not part of any Bader-atom. For the default ρ cut of 1 × 10 −3 e/Å 3 only about half of an electron is not represented by a Bader-atom. For the cutoff value needed for the calculation of the contact area, 5 × 10 −2 e/Å 3 , this number is increased to nearly 30 electrons. Although this value seems to be very large, one has to consider the total system size, which includes 3776 electrons. Thus, the relative number of "missing" electrons is below 0.8%. We also evaluated the Bader-radii R B of a single tungsten atom and a single carbon atom in a box. The value for tungsten of 2.9Å obtained for ρ cut = 5 × 10 −2 electrons perÅ 3 is more than twice as large as the empirical atomic radius, 1.35Å [47], and about 1Å larger than the calculated atomic radius of 1.93Å [48]. Reported values for the atomic radius of carbon reach from the calculated value of 0.67Å [48], over the empirical value of 0.70Å [47] to a van der Waals radius of 1.70Å [49]. As for tungsten, also the carbon Bader-atom radius for a cutoff density of 5 × 10 −2 electrons perÅ 3 is significantly larger than all these reported values at 2.30Å. This ensures that neither the tip atoms nor the surface Bader-atoms are artificially small. It also shows that it is questionable to assume that contact between two bodies is established only after the surface atoms overlap if one uses spherical atoms and traditional radii.
It is worthwhile to compare our approach to the onset of contact with the method by Mo, Turner, and Szlufarska [10], which uses the beginning of short-range interaction as a criterion for contact. As already mentioned, the distinction between long-range and short-range forces is not trivial, but it might be approximated by disabling the long-range contributions in the correlation potential of the optB86b-vdW functional. We calculated the corresponding energies at the vdW relaxed positions and fitted the data with a Morse function [50]. The interaction strength of this short-range potential at the jump to contact is 2.0% of the total potential depth which could be classified as the "beginning of the interaction". This means that our approach, at least for the system investigated here, is in accordance with the approach of reference 10. However, an interaction strength of 1%, 5% or even 10% of the short-range binding energy could also be reasonably selected as the "beginning of the interaction", each leading to different results. This highlights the advantages of using the jump to contact as the criterion for the initial point of contact, as no further assumptions are needed proceeding in this manner. Figure 4 shows a decomposition of the contact area into contributions of the tip apex (one atom) and contributions from the second tip layer (four atoms). We   FIG. 4: Ab-initio real contact area A obtained for ρ cut = 5 × 10 −2 e/Å 3 versus the distance d r in the relaxed system between the tungsten tip and graphene/Ir(111). Blue crosses give the contribution of the tip from the apex atom and green plus signs show the contribution from the second layer (four atoms). The total contact area (red circles) is given by the sum of these two contributions and is fitted by an exponential (dashed line). The inset shows the total contact area versus the static distance d s . The solid line is an exponential function resulting from equation (2) and the linear relation between d r and d s (see figure 3a).
The dotted vertical line marks the jump to contact.
find that the second layer only contributes to the total contact area for the three closest distances but is then responsible for nearly all of the increase. The dashed line in figure 4 is an exponential fit of the form to the 7 non-zero data points (red circles) with the coefficients λ r 4.2Å −1 and ∆ r 3.0Å. The factor A ∆ = 1Å 2 is included for dimensional reasons and has not been used as a fitting parameter. Although the exponential fit is not perfect, the agreement with the data is certainly reasonable, especially considering that only two fitting parameters were used. Also the point of vanishing contact is predicted well, although only points of positive contact area were considered for the fit. As there is a linear relation between d r and d s in the region where contact is established (d r = κd s + δ = 0.23d s + 1.84; see solid black line in figure 3a), it is also possible to express the contact area A through the static distance d s , which is easier accessible in experiments through the vertical displacement of the tip support. In the relation A(d s ) = A ∆ exp [−λ s (d s − ∆ s )], the decay constant is smaller than in A(d r ) with λ s = λ r κ 1.0Å −1 , while ∆ s = (∆ r − δ)/κ 5.0Å is increased compared to ∆ r , and A ∆ = 1Å 2 is the same dimensionality factor as before. This relation is plotted in the inset of figure 4. Quite naturally only the region after the jump to contact is represented well. Figure 5 shows the geometrical shape of the nonvanishing real contact area for four distances. The different colors denote different depths, ranging from ∼ 0.3Å in figure 5a to ∼ 1.4Å in figure 5b. Initially, for larger distances, the shape is rather flat and is dominated by the threefold symmetry of the graphene layer (figure 5a and 5b). As the graphene layer gets depressed towards the iridium substrate, the contact area is beginning to show a pronounced bowl shape (figure 5c), which increases in depth for decreased distance (figure 5d). Note that for the closest distance (figure 5d, d s = 1.30Å) not only the threefold symmetry of the graphene layer is visible in the center, but the edges of the bowl have the fourfold symmetry of the second tip layer.
We can also analyze how many carbon atoms are in contact with the tip for each distance and compare the contact area predicted by our approach with the results by Mo, Turner, and Szlufarska [10]. To this end we count every surface Bader atom that touches our tip as contacting, a different approach than described in the introduction, since we have no pairwise forces at our disposal. In our case, with a cross section of the simulation cell A C ∼ 262Å 2 , and 200 carbon atoms in the graphene layer the contribution per atom to the real contact area is A at = A C /200 = 1.31Å 2 . Each carbon atom has 3 nearest neighbors in d nn = 1.42Å, 6 next nearest neighbors at 2.46Å, 3 third nearest neighbors at 2.84Å, and 6 fourth nearest neighbors at 3.76Å distance. Already directly after the jump to contact at d s = 3.53Å, more than one carbon atom is in contact with the tip, although the majority of the contact is formed by the central carbon atom which is responsible for 5.90Å 2 of the total 6.82Å 2 . This is about 30% more than the contact area predicted in reference 10, with 4A at = 5.25Å 2 . The area contributed by the central carbon atom alone exceeds the value of 4A at by ∼ 12%. The next nearest and third nearest neighbors begin to play a role at d s = 2.18Å, contributing to about 7% of the total area of 16.89Å 2 . Here the method by Mo, Turner, and Szlufarska gives a very comparable area of 13A at = 17.05Å 2 . However, the 9 outermost atoms that contribute 7% to the contact area in our approach are responsible for nearly 70% of the contact area in reference 10 considering all contacting atoms equally. The situation at d s = 1.85Å is visualized in figure 6, with the method of reference 10 still giving an area of 13A at = 17.05Å 2 , while our approach yields 21.13Å 2 . Only for the two closest positions at d s = 1.51Å and d s = 1.30Å more than 13 atoms are in contact, according to the Bader partitioning, and the central 13 are still responsible for 98% and 87% of the contact area, respectively. Including the 6 fourth nearest neighbors into the model by Mo, Turner, and Szlu-farska [10], leads to 19A at = 24.92Å 2 for both of this distances, while our approach gives A(1.51) = 28.89Å 2 and A(1.30) = 35.03Å 2 . Thus, the results are comparable, but the outermost atoms are again over represented compared to our approach. Our model offers higher resolution of the real contact area and allows for a distance dependent contribution of each atom. It is important to note that our contact areas are curved and have a more or less pronounced bowl shape while Mo, Turner, and Szlufarska consider flat contact areas (see figure 6). Overall both methods show fair agreement.
Our chosen system, which has been proven to accurately model the interaction between a tungsten tip and moiré graphene on Ir(111) [17], limits our investigation to the attractive region (d r ≥ 2.24Å) and small positive loads (2.18Å ≤ d r < 2.24Å). For d r < 2.18Å, the tip forms bonds with the iridium substrate leading again to negative values of the load. Thus, it is difficult to predict the behavior of the real contact area A dependent on the load L. However, we can assume that the interaction potential E(d r ) can be approximated by a Morse potential in the vicinity of the minimum [50], where E 0 = 2.33 eV is the depth of the potential at the equilibrium position d 0 = 2.24Å, which we can get directly from our data. Thus only γ has to be fitted, resulting in γ = 4.11Å −1 . We can now derive the load Solving this for d r produces where the dimensionless function ξ(L) is Now it is possible to express the real contact area dependent on load using equations (2) and (5), which provides a power law, As γ = 4.11Å −1 and λ r = 4.19Å −1 the exponent is very close to 1, thus we arrive at a linear dependence of A on ξ(L), namely A(L) = Cξ(L) with C = 24.15Å 2 , which corresponds to an increase of A with L to the power of While we believe that our ab-initio approach using the QTAIM for calculating the real contact area is intuitive and accurate, its limitations have to be discussed as well. As the charge density is required to perform the Bader partitioning our approach is limited to system sizes where ab-initio calculations are still feasible. This limits us to a single asperity case at the moment, where only a handful of atoms interact with the surface. However, since our system was used to explain some experimental results [17,18], we are confident that our scheme is applicable to real systems, albeit only for sharp tips and low loads. With the ever increasing power of modern computers and better scaling codes, on the other hand, it might soon be feasible to calculate systems with thousands or millions of atoms and hence to study interesting phenomena such as multi-asperity contacts. Thus, our approach might be used in MD calculations as well to directly investigate the influence of the real contact area on frictional forces [51]. CONCLUSION We propose a new ab-initio approach for calculating the real contact area between a tip and a surface. We apply Bader's quantum theory of atoms in molecules (QTAIM) to determine the volumes and shapes of the atoms in contact together with their contact areas at each given distance. We define a specific density cutoff ρ cut for this partitioning to confine the Bader volumes to realistic values. This cutoff density is obtained by using the discontinuity in the d r versus d s curve to define the initial point of contact in a perspicuous way, which in the examined system occurs due to a jump to contact, commonly observed in AFM experiments. This defines a lower bound for the cutoff density which is then the optimal value, since ρ cut needs to be minimized to include the maximum number of electrons. Thus, our approach remains essentially ab-initio, as the only parameter needed can be determined from properties of the system. We believe that the jump to contact is a less ambiguous way to define the onset of contact than using a partitioning of the interaction in long-and short-range regions [10], or equating contact with repulsive interactions [11,14].
For decreasing the real tip-sample distance d r an exponential increase of the real contact area A is found. This is a combined effect of the jump to contact and the preferred distance of the tip apex relative to the surface atom below it, which first jumps up to meet the tip and then gets pressed below its equilibrium position for closer separations. As d s is linear dependent on d r , we can also express the exponential relation A(d r ) through d s , which can be better controlled in experiments than d r .

ACKNOWLEDGMENTS
The authors would like to thank A. Garhofer for providing some of the relaxed structures and helping with the selection of computational parameters, as well as F. Mittendorfer for fruitful discussions. This work was funded by the "Austrian COMET-Program" (project XTribology, no. 824187) via the Austrian Research Promotion Agency (FFG) and the Province of Niederösterreich, Vorarlberg and Wien and has been carried out within the "Excellence Centre of Tribology" (AC2T research GmbH) and at Vienna University of Technology. P.M., J.R., and G.F. acknowledge the support by the Austrian Science Fund (FWF) [SFB ViCoM F4109-N13]. This work was supported in part also by COST Action MP1303. The authors also appreciate the ample support of computer resources by the Vienna Scientific Cluster (VSC). Fig. 1 in this paper was created with the help of the VESTA code [52].