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

M Wolloch<EMAIL_ADDRESS>Institute of Applied Physics, Vienna
University of Technology, Gußhausstraße 25–25a, 1040 Vienna, Austria  Austrian
Center of Competence for Tribology, Viktor-Kaplan-Straße 2, 2700 Wiener
Neustadt, Austria  G Feldbauer  Institute of Applied Physics, Vienna
University of Technology, Gußhausstraße 25–25a, 1040 Vienna, Austria  Austrian
Center of Competence for Tribology, Viktor-Kaplan-Straße 2, 2700 Wiener
Neustadt, Austria  P Mohn  Institute of Applied Physics, Vienna University of
Technology, Gußhausstraße 25–25a, 1040 Vienna, Austria  J Redinger  Institute
of Applied Physics, Vienna University of Technology, Gußhausstraße 25–25a,
1040 Vienna, Austria  A Vernes  Institute of Applied Physics, Vienna
University of Technology, Gußhausstraße 25–25a, 1040 Vienna, Austria  Austrian
Center of Competence for Tribology, Viktor-Kaplan-Straße 2, 2700 Wiener
Neustadt, Austria

###### Abstract

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.

###### pacs:

71.15.Mb, 68.37.Ef,68.37.Ps,73.63.Rt,68.55.ap

##  I 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 Bowden and Tabor
(1939). 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
Hertz (1881). 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 Johnson  _et al._ (1971). 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 Derjaguin  _et al._ (1975). 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 Tabor (1977). Additional work on this unification
was done by Maugis using Dugdale potentials Maugis (1992). 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 Carpick  _et al._ (1999). 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)\sim 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 Luan and Robbins
(2005). 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 Dong  _et al._ (2013).
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 Luan and Robbins (2005); Mo  _et al._
(2009); Cheng and Robbins (2010); Eder  _et al._ (2011, 2014). 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 Burnham  _et al._
(1991); Cheng and Robbins (2010). 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 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 Mo  _et al._ (2009). For this large scale finite
temperature (300 K) MD study they used a reactive empirical bond-order (REBO)
potential Brenner  _et al._ (2002) 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
(Mo and Szlufarska, 2010). Now an atomic contact area $A_{\text{at}}$ is
attributed to every chemically interacting atom of the tip leading to the
total real contact area $A_{\text{real}}=N_{\text{at}}A_{\text{at}}$, with
$N_{\text{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_{\text{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_{\text{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 Voloshina  _et al._ (2013); Garhofer
(2013). 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.

##  II Computational methods

The simulation cell consists of 4 layers of a $9\times 9$ iridium(111)
substrate covered by a $10\times 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 Klimeš  _et al._ (2010,
2011), $a_{\mathrm{Ir}}=2.735\,\mathrm{\AA}$ and
$a_{\mathrm{Gr}}=2.465\,\mathrm{\AA}$, are close to the experimental values of
2.71 Å and 2.46 Å, respectively. The mismatch of the structure is very small
at $10a_{\mathrm{Gr}}-9a_{\mathrm{Ir}}=0.015\,\mathrm{\AA}$ Garhofer (2013).
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 ($\sim$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.

(a) Side view

(b) Top view

Figure 1: Side (a) and top (b) view of a tungsten tip on graphene/Ir(111).
Iridium atoms are shown in yellow, carbon in brown and tungsten in grey.

All calculations were performed within density functional theory (DFT)
employing the Vienna Ab-Initio Simulation Package VASP Kresse and Hafner
(1993, 1994); Kresse and Furthmüller (1996a, b) using the Projector Augmented-
Wave (PAW) method Blöchl (1994); Kresse and Joubert (1999). To include van der
Waals (vdW) forces, which are relevant in this system, the optB86b-vdW
functional was employed Klimeš  _et al._ (2010, 2011). This vdW density
functional has been applied to a wide range of materials and proven to be of
good accuracy Chakarova-Käck  _et al._ (2006); Sony  _et al._ (2007); Carrasco
_et al._ (2011); Mittendorfer  _et al._ (2011); Graziano  _et al._ (2012);
Bedolla  _et al._ (2014a, b). The Brillouin zone sampling was performed on a
$\Gamma$-centered $3\times 3\times 1$ k-grid, with a smearing of 0.1 eV using
the method of Methfessel and Paxton to first order Methfessel and Paxton
(1989). To ensure good accuracy for the Bader-partitioning scheme, the
electronic charge density was calculated on a dense mesh of $432\times
432\times 448$ points in the simulation cell. Electronic energies were
converged to $10^{-6}\,\mathrm{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 Garhofer (2013), 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 $\rho$ in our simulation cell into
single atoms we use Bader’s quantum theory of atoms in molecules (QTAIM) Bader
(1990). In contrast to other similar approaches Hirshfeld (1977); Bultinck
_et al._ (2007); Bickelhaupt  _et al._ (1996); Fonseca Guerra  _et al._
(2004); Mayer and Salvador (2004); Becke (1988), 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 $\rho$ and is given as Bader (1985),

$\nabla\rho({\bf r})\cdot{\bf n}({\bf r})=0\quad\forall{\bf r}\in S({\bf
r})\quad.$ (1)

Here $S$ is the boundary surface of the atom and ${\bf n}(\bf r)$ is the unit
vector normal to this surface. The condition states that the flux of the
gradient field of the charge density, $\nabla\rho$, 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 Henkelman  _et al._ (2006); Sanville  _et al._ (2007); Tang  _et al._
(2009).

##  III 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 $\rho$ 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 $\rho_{\mathrm{cut}}$ cannot be chosen
arbitrarily, since it directly influences the contact area.

A possibility to extract a value for $\rho_{\mathrm{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 Mo  _et al._ (2009). This
procedure defines a density cutoff $\rho_{\mathrm{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.

Figure 2: Sketch of possible tip-surface interactions. If the tip is far away
from the surface (middle panel), 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.

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\,\mathrm{\AA}$ and $d_{s}=3.53\,\mathrm{\AA}$ (see figure 3(a)),
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 3(b)). 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}\simeq 3.6\,\mathrm{\AA}$
and to tune the value of the density cutoff $\rho_{\mathrm{cut}}$ accordingly.
Note that an infinitely stiff cantilever is assumed in our calculations, as
the uppermost atoms of our tip are kept 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.

(a)

(b)

Figure 3: (a) Distance $d_{r}$ in the relaxed system (red crosses) versus the
(static) distance $d_{s}$ between a rigid tungsten tip approaching
graphene/Ir(111). The dashed line gives $d_{r}=d_{s}$. Attractive interactions
cause a jump to contact between $d_{s}=3.65\,\mathrm{\AA}$ and
$d_{s}=3.53\,\mathrm{\AA}$, which is marked with a vertical dotted line. After
the jump to contact the relation between $d_{r}$ and $d_{s}$ is also linear
(solid black line). This line crosses the $d_{r}=d_{s}$ line at the
equilibrium point, where the graphene layer has relaxed back into its original
shape. A hypothetical curve for a purely repulsive interaction is sketched by
the dashed-dotted line to illustrate that the method is also viable if no jump
to contact is occurring in the system. (b) Displacement $D_{C}$ of the carbon
atom situated directly below the tip apex with respect to its initial
position, versus the distance $d_{s}$. The jump to contact is again marked by
a vertical dotted line.

Once $\rho_{\mathrm{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
Henkelman  _et al._ (2006); Sanville  _et al._ (2007); Tang  _et al._ (2009),
up to a cutoff of 10-1 e/Å3. While the default value of
$\rho_{\mathrm{cut}}=10^{-3}\,\mathrm{e/\AA^{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\,\mathrm{\AA}$, for
a value of $\rho_{\mathrm{cut}}\sim 10^{-2}\,\mathrm{e/\AA^{3}}$. Obviously, a
very high $\rho_{\mathrm{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
$\rho_{\mathrm{cut}}$ ranging from $7.5\times 10^{-2}$ e/Å3 to $1.0\times
10^{-2}$ e/Å3 in order to find the lowest value that still satisfies these
conditions, resulting in an optimal value of $\rho_{\mathrm{cut}}=5\times
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 $\rho_{\mathrm{cut}}$. However, the obtained cutoff density of
$\rho_{\mathrm{cut}}=5\times 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 111Calculations where carried out on the Vienna Scientific Cluster,
typically using 16 computer nodes interconnected by an QDR InfiniBand network,
each equipped with two 8-core AMD Opteron 6132HE 2.2 GHz CPUs and 32 GB of ECC
DDR3 RAM. Computation times on this cluster, including pre-relaxation with a
GGA functional, full relaxations with the optB86b vdW functional and final
static calculations with the finer mesh to generate proper charge densities
for the Bader partitioning, took up to 100 hours or 25600 corehours for each
tip position., 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\times 10^{-2}$ electrons
per Å3, which is approximately the same as for the graphene/Ir(111) system.

As the optimized cutoff of $5\times 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 $\rho_{\mathrm{cut}}$ of $1\times
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\times 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_{\mathrm{B}}$ of a single tungsten atom and a single carbon atom in a box.
The value for tungsten of 2.9 Å obtained for $\rho_{\mathrm{cut}}=5\times
10^{-2}$ electrons per Å3 is more than twice as large as the empirical atomic
radius, 1.35 Å Slater (1964), and about 1 Å larger than the calculated atomic
radius of 1.93 Å Clementi  _et al._ (1967). Reported values for the atomic
radius of carbon reach from the calculated value of 0.67 Å Clementi  _et al._
(1967), over the empirical value of 0.70 Å Slater (1964) to a van der Waals
radius of 1.70 Å Bondi (1964). As for tungsten, also the carbon Bader-atom
radius for a cutoff density of $5\times 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 Mo  _et al._ (2009), 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 Morse (1929). 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: Ab-initio real contact area $A$ obtained for
$\rho_{\mathrm{cut}}=5\times 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 3(a)). The dotted vertical line marks the jump to contact.

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 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

$A(d_{r})=A_{\Delta}e^{-\lambda_{r}\left(d_{r}-\Delta_{r}\right)}\quad,$ (2)

to the 7 non-zero data points (red circles) with the coefficients
$\lambda_{r}\simeq 4.2\,\mathrm{\AA^{-1}}$ and $\Delta_{r}\simeq
3.0\,\mathrm{\AA}$. The factor $A_{\Delta}=1\,\mathrm{\AA^{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}=\kappa d_{s}+\delta=0.23d_{s}+1.84$; see solid black line
in figure 3(a)), 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_{\Delta}\exp\left[-\lambda_{s}\left(d_{s}-\Delta_{s}\right)\right]$,
the decay constant is smaller than in $A(d_{r})$ with
$\lambda_{s}=\lambda_{r}\kappa\simeq 1.0\,\mathrm{\AA^{-1}}$, while
$\Delta_{s}=(\Delta_{r}-\delta)/\kappa\simeq 5.0\,\mathrm{\AA}$ is increased
compared to $\Delta_{r}$, and $A_{\Delta}=1\,\mathrm{\AA^{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 non-vanishing real contact area
for four distances. The different colors denote different depths, ranging from
$\sim 0.3\,\mathrm{\AA}$ in figure 5(a) to $\sim 1.4\,\mathrm{\AA}$ in figure
5(b). Initially, for larger distances, the shape is rather flat and is
dominated by the threefold symmetry of the graphene layer (figure 5(a) and
5(b)). As the graphene layer gets depressed towards the iridium substrate, the
contact area is beginning to show a pronounced bowl shape (figure 5(c)), which
increases in depth for decreased distance (figure 5(d)). Note that for the
closest distance (figure 5(d), $d_{s}=1.30\,\mathrm{\AA}$) 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.

(a) $d_{s}=3.53\,\mathrm{\AA}$

(b) $d_{s}=3.14\,\mathrm{\AA}$

(c) $d_{s}=1.85\,\mathrm{\AA}$

(d) $d_{s}=1.30\,\mathrm{\AA}$

Figure 5: Contact formed by lowering a tungsten tip onto a graphene/Ir(111)
surface. The contact area increases from (a) 6.8 Å2, over (b) 9.1 Å2, and (c)
21.13 Å2, to (d) 35.0 Å2. Different depths of the curved contact areas are
coded by color contours with the lowest value set to zero. Please note the
different color bars and axes scaling in each panel.

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 Mo  _et al._ (2009). 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_{\mathrm{C}}\sim 262\,\mathrm{\AA^{2}}$, and 200 carbon atoms in the
graphene layer the contribution per atom to the real contact area is
$A_{\mathrm{at}}=A_{\mathrm{C}}/200=1.31\,\mathrm{\AA^{2}}$. Each carbon atom
has 3 nearest neighbors in $d_{\mathrm{nn}}=1.42\,\mathrm{\AA}$, 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\,\mathrm{\AA}$, 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_{\mathrm{at}}=5.25\,\mathrm{\AA^{2}}$. The area contributed by the central
carbon atom alone exceeds the value of $4A_{\mathrm{at}}$ by $\sim 12\%$. The
next nearest and third nearest neighbors begin to play a role at
$d_{s}=2.18\,\mathrm{\AA}$, 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_{\mathrm{at}}=17.05\,\mathrm{\AA^{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\,\mathrm{\AA}$ is
visualized in figure 6, with the method of reference 10 still giving an area
of $13A_{\mathrm{at}}=17.05\,\mathrm{\AA^{2}}$, while our approach yields
21.13 Å2. Only for the two closest positions at $d_{s}=1.51\,\mathrm{\AA}$ and
$d_{s}=1.30\,\mathrm{\AA}$ 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 Szlufarska Mo  _et al._ (2009), leads to
$19A_{\mathrm{at}}=24.92\,\mathrm{\AA^{2}}$ for both of this distances, while
our approach gives $A(1.51)=28.89\,\mathrm{\AA^{2}}$ and
$A(1.30)=35.03\,\mathrm{\AA^{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.

Figure 6: Top view (left) and side view (right) of the real contact area
resulting from our ab-initio approach using Bader atoms (color code depending
on height) for a distance of $d_{s}=1.85\,\mathrm{\AA}$ (see figure 5(c)),
compared to the flat contact area from the model by Mo, Turner, and Szlufarska
(green) Mo  _et al._ (2009), Carbon and tungsten atoms are sketched as red and
black dots, respectively.

Our chosen system, which has been proven to accurately model the interaction
between a tungsten tip and moiré graphene on Ir(111) Voloshina  _et al._
(2013), limits our investigation to the attractive region ($d_{r}\geq
2.24\,\mathrm{\AA}$) and small positive loads ($2.18\,\mathrm{\AA}\leq
d_{r}<2.24\,\mathrm{\AA}$). For $d_{r}<2.18\,\mathrm{\AA}$, 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 Morse (1929),

$E_{\mathrm{M}}\left(d_{r}\right)=E_{0}\left\\{\left[1-e^{-\gamma\left(d_{r}-d_{0}\right)}\right]^{2}-1\right\\}\quad,$
(3)

where $E_{0}=2.33\,\mathrm{eV}$ is the depth of the potential at the
equilibrium position $d_{0}=2.24\,\mathrm{\AA}$, which we can get directly
from our data. Thus only $\gamma$ has to be fitted, resulting in
$\gamma=4.11\,\mathrm{\AA^{-1}}$. We can now derive the load $L=-\partial
E_{\mathrm{M}}/\partial d_{r}$ yielding

$L\left(d_{r}\right)=-2\gamma
E_{0}\left[1-e^{-\gamma\left(d_{r}-d_{0}\right)}\right]e^{-\gamma\left(d_{r}-d_{0}\right)}\quad.$
(4)

Solving this for $d_{r}$ produces

$d_{r}=d_{0}-\dfrac{1}{\gamma}\ln\xi\left(L\right)\quad,$ (5)

where the dimensionless function $\xi(L)$ is

$\xi\left(L\right)=\dfrac{1\pm\sqrt{1-4u}}{2}~{},\qquad\mathrm{with}\qquad
u=\dfrac{L}{2\gamma E_{0}}\quad.$ (6)

Now it is possible to express the real contact area dependent on load using
equations (2) and (5), which provides a power law,

$A\left(L\right)=A_{0}e^{-\lambda_{r}\left(d_{0}-\Delta_{r}\right)}\,\left[\xi\left(L\right)\right]^{\frac{\lambda_{r}}{\gamma}}\quad.$
(7)

As $\gamma=4.11\,\mathrm{\AA^{-1}}$ and $\lambda_{r}=4.19\,\mathrm{\AA^{-1}}$
the exponent is very close to 1, thus we arrive at a linear dependence of $A$
on $\xi(L)$, namely $A(L)=C\xi(L)$ with $C=24.15\,\mathrm{\AA^{2}}$, which
corresponds to an increase of $A$ with $L$ to the power of $\frac{1}{2}$.

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 Voloshina  _et
al._ (2013); Garhofer (2013), 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 Vernes  _et al._ (2012).

##  IV 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 $\rho_{\mathrm{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
$\rho_{\mathrm{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 Mo  _et al._ (2009), or equating contact with repulsive interactions
Burnham  _et al._ (1991); Cheng and Robbins (2010).

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 Momma and Izumi (2011).

## References

  * Bowden and Tabor (1939) F. P. Bowden and D. Tabor, Royal Society of London Proceedings Series A 169, 391 (1939).
  * Hertz (1881) H. Hertz, Journal für die reine und angewandte Mathematik 92, 156 (1881).
  * Johnson  _et al._ (1971) K. L. Johnson, K. Kendall, and A. D. Roberts, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 324, 301 (1971) .
  * Derjaguin  _et al._ (1975) B. Derjaguin, V. Muller, and Y. Toporov, Journal of Colloid and Interface Science 53, 314 (1975).
  * Tabor (1977) D. Tabor, Journal of Colloid and Interface Science 58, 2 (1977), international Conference on Colloids and Surfaces.
  * Maugis (1992) D. Maugis, Journal of Colloid and Interface Science 150, 243 (1992).
  * Carpick  _et al._ (1999) R. W. Carpick, D. Ogletree, and M. Salmeron, Journal of Colloid and Interface Science 211, 395 (1999).
  * Luan and Robbins (2005) B. Luan and M. O. Robbins, Nature 435, 929 (2005).
  * Dong  _et al._ (2013) Y. Dong, Q. Li, and A. Martini, Journal of Vacuum Science and Technology A: Vacuum, Surfaces, and Films 31, 030801 (2013).
  * Mo  _et al._ (2009) Y. Mo, K. T. Turner, and I. Szlufarska, Nature 457, 1116 (2009).
  * Cheng and Robbins (2010) S. Cheng and M. Robbins, Tribology Letters 39, 329 (2010).
  * Eder  _et al._ (2011) S. Eder, A. Vernes, G. Vorlaufer, and G. Betz, Journal of Physics: Condensed Matter 23, 175004 (2011).
  * Eder  _et al._ (2014) S. Eder, A. Vernes, and G. Betz, Comput. Phys. Commun. 185, 217 (2014).
  * Burnham  _et al._ (1991) N. A. Burnham, R. J. Colton, and H. M. Pollock, Journal of Vacuum Science & Technology A 9, 2548 (1991).
  * Brenner  _et al._ (2002) D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, Journal of Physics Condensed Matter 14, 783 (2002).
  * Mo and Szlufarska (2010) Y. Mo and I. Szlufarska, Phys. Rev. B 81, 035405 (2010).
  * Voloshina  _et al._ (2013) E. N. Voloshina, E. Fertitta, A. Garhofer, F. Mittendorfer, M. Fonin, A. Thissen, and Y. S. Dedkov, Sci. Rep. 3, 1072 (2013).
  * Garhofer (2013) A. Garhofer,  _Ab initio studies of graphene-metal interfaces_ , Ph.D. thesis, Vienna University of Technology (2013).
  * Klimeš  _et al._ (2010) J. Klimeš, D. R. Bowler, and A. Michaelides, Journal of Physics: Condensed Matter 22, 022201 (2010).
  * Klimeš  _et al._ (2011) J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Comput. Mat. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Chakarova-Käck  _et al._ (2006) S. Chakarova-Käck, E. Schröder, B. Lundqvist, and D. Langreth, Physical Review Letters 96, 146107 (2006).
  * Sony  _et al._ (2007) P. Sony, P. Puschnig, D. Nabok, and C. Ambrosch-Draxl, Phys. Rev. Lett. 99, 176401 (2007).
  * Carrasco  _et al._ (2011) J. Carrasco, B. Santra, J. Klimeš, and A. Michaelides, Physical Review Letters 106, 026101 (2011).
  * Mittendorfer  _et al._ (2011) F. Mittendorfer, A. Garhofer, J. Redinger, J. Klimeš, J. Harl, and G. Kresse, Physical Review B 84, 201401(R) (2011).
  * Graziano  _et al._ (2012) G. Graziano, J. Klimeš, F. Fernandez-Alonso, and A. Michaelides, Journal of Physics: Condensed Matter 24, 424216 (2012).
  * Bedolla  _et al._ (2014a) P. O. Bedolla, G. Feldbauer, M. Wolloch, S. J. Eder, N. Dörr, P. Mohn, J. Redinger, and A. Vernes, The Journal of Physical Chemistry C 118, 17608 (2014a) .
  * Bedolla  _et al._ (2014b) P. O. Bedolla, G. Feldbauer, M. Wolloch, C. Gruber, S. J. Eder, N. Dörr, P. Mohn, J. Redinger, and A. Vernes, The Journal of Physical Chemistry C 118, 21428 (2014b) .
  * Methfessel and Paxton (1989) M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989).
  * Bader (1990) R. F. W. Bader,  _Atoms in Molecules - A Quantum Theory_ (University of Oxford Press, Oxford, 1990).
  * Hirshfeld (1977) F. Hirshfeld, Theoretica chimica acta 44, 129 (1977).
  * Bultinck  _et al._ (2007) P. Bultinck, C. Van Alsenoy, P. W. Ayers, and R. Carbó-Dorca, The Journal of Chemical Physics 126, 144111 (2007).
  * Bickelhaupt  _et al._ (1996) F. M. Bickelhaupt, N. J. R. van Eikema Hommes, C. Fonseca Guerra, and E. J. Baerends, Organometallics 15, 2923 (1996) .
  * Fonseca Guerra  _et al._ (2004) C. Fonseca Guerra, J.-W. Handgraaf, E. J. Baerends, and F. M. Bickelhaupt, Journal of Computational Chemistry 25, 189 (2004).
  * Mayer and Salvador (2004) I. Mayer and P. Salvador, Chemical Physics Letters 383, 368 (2004).
  * Becke (1988) A. D. Becke, The Journal of Chemical Physics 88, 2547 (1988).
  * Bader (1985) R. F. W. Bader, Accounts of Chemical Research 18, 9 (1985) .
  * Henkelman  _et al._ (2006) G. Henkelman, A. Arnaldsson, and H. Jónsson, Comput. Mater. Sci. 36, 354 (2006).
  * Sanville  _et al._ (2007) E. Sanville, S. D. Kenny, R. Smith, and G. Henkelman, Journal of Computational Chemistry 28, 899 (2007).
  * Tang  _et al._ (2009) W. Tang, E. Sanville, and G. Henkelman, Journal of Physics: Condensed Matter 21, 084204 (2009).
  * Note (1) Calculations where carried out on the Vienna Scientific Cluster, typically using 16 computer nodes interconnected by an QDR InfiniBand network, each equipped with two 8-core AMD Opteron 6132HE 2.2 GHz CPUs and 32 GB of ECC DDR3 RAM. Computation times on this cluster, including pre-relaxation with a GGA functional, full relaxations with the optB86b vdW functional and final static calculations with the finer mesh to generate proper charge densities for the Bader partitioning, took up to 100 hours or 25600 corehours for each tip position.
  * Slater (1964) J. C. Slater, Journal of Chemical Physics 41, 3199 (1964).
  * Clementi  _et al._ (1967) E. Clementi, D. L. Raimondi, and W. P. Reinhardt, The Journal of Chemical Physics 47, 1300 (1967).
  * Bondi (1964) A. Bondi, The Journal of Physical Chemistry 68, 441 (1964) .
  * Morse (1929) P. M. Morse, Phys. Rev. 34, 57 (1929).
  * Vernes  _et al._ (2012) A. Vernes, S. Eder, G. Vorlaufer, and G. Betz, Faraday Discuss. 156, 173 (2012).
  * Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).