Abstract
In this work, singleasperity contact mechanics is investigated for positive and negative work of adhesion Δγ. In the latter case, finiterange repulsion acts in addition to hardwall constraints. This constitutes a continuum model for a contact immersed in a strongly wetting fluid, which can only be squeezed out in the center of the contact through a sufficiently large normal load F_{N}. As for positive work of adhesion, two stable solutions can coexist in a finite range of normal loads. The competing solutions can be readily interpreted as contacts with either a loadbearing or a squeezedout fluid. The possibility for coexistence and the subsequent discontinuous wetting and squeezeout instabilities depend not only on the Tabor coefficient μ_{T} but also on the functional form of the finiterange repulsion. For example, coexistence and discontinuous wetting or squeezeout do not occur when the repulsion decreases exponentially with distance. For positive work of adhesion, the normal displacement mainly depends on F_{N}, Δγ, and μ_{T} but – unlike the contact area – barely on the functional form of the finiterange attraction. The results can benefit the interpretation of atomic force microscopy in liquid environments and the modeling of multiasperity contacts.
Introduction
The continuum description of singleasperity contact mechanics has received much attention in the last few decades. This is in large parts because it describes forcedisplacement curves rather accurately down to nanometer scales relevant to atomic force microscopy (AFM) [13]. The contributions to the linear elasticity of (frictionless) singleasperity contacts most central to this work are the following: Hertz [4] solved the contact mechanics of a parabolic tip pressed against a substrate for hardwall repulsion. He found that the contact area A_{c} and the separation between the two solids d both disappear continuously with as the normal load squeezing the two solids together, F_{N}, approaches zero. Derjaguin, Muller, and Toporov (DMT) [5] included adhesion as a longrange interaction, in which case adhesion effectively acts as a normal load that is independent of the contactradius a_{c}. Johnson, Kendall, and Roberts (JKR) [6] solved the problem for shortrange adhesion, for which the attractive surface forces act directly at the contact line. Unlike Hertz and DMT, JKR found finite values for A_{c} and d at pull off. Tabor [7] introduced a dimensionless parameter, μ_{T}, now known as Tabor coefficient, allowing one to estimate if a contact is closer to DMT or to JKR theory. He actually recognized that DMT and JKR describe the opposite limits of long and shortrange forces, respectively. This had not been known before but was soon confirmed in numerical simulations by Muller, Yushenko, and Derjaguin [8]. Lastly, Maugis [9] used a cohesivezone model introduced by Dugdale (MD) and found analytical solutions for intermediaterange adhesion at arbitrary values of μ_{T}.
Although singleasperity, linearlyelastic, adhesive contacts mechanics is a rather mature field [10], two key issues remain worth addressing: First, only few studies have considered the case of negative work of adhesion [11,12], Δγ < 0, specifically finiterange repulsion between two surfaces acting in addition to hardwall repulsion. In particular, the concept of the Tabor coefficient has not yet been applied to negative Δγ. Therefore, I investigate if there are different regimes for Δγ < 0 as is the case for contacts with Δγ > 0, which are classified as JKR for large μ_{T} and as DMT for small μ_{T}. This includes a study of which parameters determine the behavior near squeezeout as well as a comparison to the behavior near pulloff for Δγ > 0. For the latter, it is straightforward to deduce from established results how the a_{c}(F_{N}) relation depends on the Tabor coefficient in the DMT and the JKR limit. Specifically, a_{c} − a_{p} (F_{N} + F_{p})^{κ} for F_{N} ≥ −F_{p}, where F_{p} and a_{p} are pulloff force and pulloff radius, respectively. They both depend on μ_{T} just like the exponent κ, e.g., a_{p} > 0 and κ = 1/2 in the JKR limit, while a_{p} = 0 and κ = 1/3 for DMT. In the present comparison of squeezeout (finiterange repulsion) versus pulloff (finiterange attraction), I also study whether the exponent κ changes continuously between JKR and DMT or discontinuously – as assumed implicitly in the Carpick–Ogletree–Salmeron (COS) model [1].
The second motivation for this paper is that it has not yet been investigated sufficiently how the (precise) functional form for adhesive interactions affects contact mechanics – assuming that all continuum parameters, from normal load to Tabor coefficient, are identical. It is only established that there is little sensitivity in the limits of large and zero Tabor coefficients. Yet, when studying contactmechanics between macroscopic, adhesive, rough surfaces in the context of continuummechanics, one would want to know how to best reach the JKR limit, which appears to be the relevant limit for that application. For example, it is used implicitly in Persson theory for nominally flat, adhesive surfaces [13]. In fact, the current work was initiated by the desire to add adhesion into a Green’s function molecular dynamics (GFMD) code used to model the contact between rough surfaces. To model adhesion, one needs to identify a functional form for the finiterange surface forces allowing one to reach the JKR limit in an efficient and wellcontrolled fashion. It quickly became clear that doing a clean job is not entirely trivial and that modeling singleasperity contacts ought to be better understood first and moreover is interesting in its own right.
The remainder of this paper is organized as follows: I first introduce the model, sketch the numerical methods and discuss difficulties arising in simulations in the limit of large and small Tabor coefficients. Next, I present a brief dimensional analysis motivating the commonly used unit system and the Tabor coefficient. The result section opens with adhesive contacts. There, I reproduce some established results and investigate how sensitive results are on the details of the interaction model. That section also contains a comparison to and an asymptotic analysis of the MD model motivating some minor modifications of the empirical COS equations [1]. Next, results on positive adhesion are presented before the major results are summarized and conclusions are drawn.
Results and Discussion
Definition of the model
In this section, the singleasperity contact problem is introduced. As shown in Figure 1, we consider a stiff, ideallyflat wall positioned in the xy plane at z = 0 and a linearlyelastic tip of parabolic shape. Its undeformed surface is given by
where R is the radius of curvature and the inplane distance of the center of the tip from the origin of the coordinate system. The elastic displacement of the tip, u(x,y) is formally a function of both inplane coordinates, although the equilibrium solutions only depend on r. The gap g(x,y) indicates the distance between the deformed tip and the substrate, i.e.,
It is furthermore assumed that the tip cannot penetrate the substrate. This can be stated as a nonholonomic boundary condition
Alternatively, one can formally introduce a shortrange, hardwall repulsion [14]
where f_{r} is an arbitrary positive constant of unit force per area. Note that the integrand on the r.h.s. of Equation 4 is zero for finite gaps while it diverges for negative gaps. Depending on the problem, it can be more convenient to use either the nonholonomic boundary condition or the energybased description of the shortrange repulsion.
This work also considers finiterange adhesive or finiterange repulsive interactions V_{fr}, which only depend on the local gap. The default expression for it is:
where Δγ is the work of adhesion per surface area that is obtained when a flat tip touches the substrate in a continuum description. The choice of the functional dependence of V_{fr} is not motivated by the true functional form for the interactions between any two real solids, but for the moment being, it is a matter of convenience. Alternative interaction models for the integrand are introduced in a seperate section.
An important property of all models for V_{fr} is that the interaction is characterized by a prefactor representing the work of adhesion and a single length scale z_{0}. The choice of the latter allows one to localize adhesive stress near the contact line through z_{0} → 0 or to extend the adhesive interactions to radii much exceeding the contact radius a_{c} for z_{0} → ∞. Of course, z_{0} can take any value between zero and infinity so that intermediaterange interactions can be modeled as well. Note that V_{fr} and V_{sr} are qualitatively different: The prefactor of the shortrange potential is formally zero, because f_{r} is finite and z_{r} very small. In other words, z_{r} represents the size of an “infinitesimallysmall” atom whose size is irrelevant on a continuum scale. In contrast, the prefactor of the finiterange potential is considered finite as well as the range of interaction z_{0}. It represents a “collective” length scale, such as the decay length of density oscillations in the fluid [15] or the radius of gyration of a polymer.
The displacement u(x,y) and other fields (gap and stress) will be expressed not only in real space, but also in Fourier space. This is done most conveniently by using inplane periodic boundary conditions. The respective boundaries lie at x or y = ±L/2, which should be chosen such that L (the linear dimension of the simulation cell) is much greater than the linear dimension of the contact zone. The latter includes the contact and the area of nonnegligible adhesive (or finiterange repulsive forces) stresses. The following convention for the Fourier transform shall be employed
where the wave vector components satisfy q_{α} = 2πn/L, A = L^{2} is the integration domain, and n an integer number. With these definitions, one can express the elastic energy of the deformed tip (in the smallslope approximation) as
where E is the effective modulus, E = E_{Y}/(1 − ν^{2}), E_{Y} being the Young’s modulus and ν the Poisson ratio. The convention of using the symbol E^{*} for the effective modulus is abandoned for clarity, because primes will be used later to indicate scaled coordinates and scaled parameters.
Since can be interpreted as the center of mass mode, the effect of a load (or normal force) exerted on the tip leads to an energy
When solving the contact problem, one seeks to minimize the total energy
with respect to u(x,y), i.e., the solution u_{0}(r) must satisfy
for each r. Here, δ indicates a functional derivative. In a discrete representation of the problem, r is an index so that the functional derivative in Equation 11 has to be replaced by a partial derivative.
Alternative interaction models
Throughout this paper, different functional forms for the finiterange interactions between surfaces are considered in addition to the “default” or “exponential” model introduced in Equation 5. Functions similar to the ones used in this work have already been employed for the simulation of mode I fracture or delamination. Depending on the authors, the functions are called the cohesive zone model [16], the crack evolution function [17], or the tractionseparation relation [18]. However, it is not clear how the results obtained for mode I fracture geometries relate to Hertzian contacts. This is the main reason why the results obtained within the cohesive zone model cannot be compared in a straightforward fashion to those of the current study.
The additional models in the current work replace the integrand on the r.h.s. of Equation 5 with the following expressions:
where Θ(···) denotes the Heavyside step functions, which is one for positive arguments and zero otherwise. The interaction potentials and their first derivatives are shown in Figure 2.
All expressions take the same value, −Δγ, for the adhesive energy when the two surfaces touch, i.e., in each case the work of adhesion is Δγ. In this sense, all four models produce the same continuum limit. However, in two models, namely the Gauss and the van der Waals (VDW) integrands, the derivative goes to zero when the two surfaces touch, while it remains finite in the exponential model and the MD model.
As stated before, none of the models are supposed to be highly realistic representations of surface forces, although each model may have its justification. In particular the exponential model follows from the argument that longrange density correlations in fluids are governed by a single length [15]. In a highdensity fluid, the correlation length becomes complex [15], which then leads to layering transitions as discussed recently [19]. The VDW model might approximate the longrange van der Waals interactions in a way that Δγ reflects the Hamaker constant. Depending on the confined system in question, other effective interactions might be possible. However, all models represent the feature that surfaces repel upon close approach (i.e., when atoms from opposing surfaces bump into each other, which is implemented through the hardwall repulsion) and that attraction – or additional repulsion – may occur at finite distance. Continuous shortrange repulsive forces are not used here. Doing so would complicate the definition of contact and thus contact radius, which has remained controversial for systems without hardwall or harddisk interactions [2023]. Lastly, the equations to be solved would become very stiff and thus the simulation inefficient if the hardwall constraint was replaced by potentials with large curvature.
In the context of the squeezeout of fluids, the MD and the exponential model might not be physically meaningful for small ratios of g/z_{0}: when one flat solid is placed on top of another flat solid with an infinitesimally small external load (in the absence of a fluid), the two solids would repel, although they “cannot know”, away from a contact line, that a fluid wants to penetrate. Nonetheless, the exponential model has been used in early study of squeezeout of fluids [12]. Forces between two (flat) surfaces in the Gauss and the VDW model are zero either for intimate mecanical contact or at infinite separation.
Dimensional analysis
In this section, I present a simple dimensional analysis of the contact problem. The result of the analysis is a meaningful set of units, which, in similar form, has already been established by Maugis [9]. However, in the present analysis, units are not motivated from the solutions but rather straight from the beginning, i.e., by the expressions defining the model. This is why Maugis’ and the present units differ by dimensionless prefactors, which, however, always turn out close to unity. In the subsequent derivation, it is not necessary to know the precise functional dependence of the finiterange forces.
Assume we know the solution u_{0}(r) minimizing V_{tot} for a given set of parameters defining our model, i.e., u(r) solves the contact problem for a specific set of values for E, R, F_{N}, Δγ, and z_{0}. It is then possible to “recycle” the shape of the function u_{0}(r) to solve a different problem defined by different parameters E′, R′, Δγ′, and . Specifically, if each individual summand of is identical to the equivalent term in V_{tot}[u(r)] (up to a multiplicative constant, which can be set to one), then minimizes given that u_{0}(r) minimizes V_{tot}[u(r)].
The transformation, in which inplane coordinates are scaled as r′ ≡ sr and normal coordinates are scaled as z′ = tz leaves the shape of the solution unchanged. Of course, z(r) and thus g(r) must be transformed the same way as u(r). Therefore, the radius of curvature of the “new” tip is R′ = s^{2}R/t.
Let us investigate how one has to alter each individual term contributing to so that it matches its equivalent in V_{tot}[u]. (i) The hardwall repulsive energy V_{fr} is unproblematic. It disappears for the old and the new solution, because of the limit z_{r} → 0, i.e., . (ii) To recycle the V_{fr} calculation, we need to set = tz_{0}. Keeping in mind that A′ = s^{2}A, where A = L^{2} is the original integration domain, it follows that = s^{2}(Δγ′/Δγ)V[u]. (iii) For the calculation of the elastic energy, it is useful to keep in mind that q′ = q/s and that A′ = s^{2}A. This means that (The integer indices enumerating the wave vectors are identical for the original and the new domain.) (iv) Lastly, the loadrelated energy becomes In summary, we can recycle our solution with the following substitutions
Let us first consider the case of no external force, F_{N} = 0, so that we investigate the “intrinsic” system properties. If we use E as the unit of pressure, which is done until further notice, all three remaining parameters defining the system can be expressed to be of unit length, i.e, R, z_{0}, and Δγ/E. Wether a potential should be classified as short or longranged can only depend on a nondimensionalized interaction length. This means that z_{0} has to be expressed in the two remaining units of length (R and Δγ/E) such that the dimensionless ratio
remains unchanged under a scaling transformation.
Let us now chose the radius of curvature as the unit of length, or, alternatively, consider only those scaling transformations that leave R constant. This can be achieved by setting t = s^{2}, which maps an infinite parabola (x → z = x^{2}) onto itself (sx → s^{2}z). I note in passing that such a transformation might not be meaningful for a scaling analysis of the contact mechanics of randomly rough surfaces, which will be presented elsewhere.
Inserting the relevant equalities from Equation 13–Equation 17 reveals that choosing α = 2/3 leaves the expression in Equation 18 constant. As a consequence, the range of influence of the adhesive term is best quantified by the term
where μ_{T} is known as Tabor coefficient – up to a dimensionless, multiplicative constant. It remains invariant under all scaling transformation in Equation 13–Equation 17 leaving the radius of curvature unchanged.
From Equation 19 one can see that we need to send z_{0} → w^{2/3}z_{0} in order to keep the Tabor coefficient constant when changing Δγ/E at constant R to wΔγ/E. This in turn implies a transformation of x → w^{1/3}x for the inplane coordinates, because R is supposed to remain unchanged. It follows that a_{c}(F_{N} = 0) → w^{1/3}a_{c}(F_{N} = 0) and thus a_{c} (Δγ/E)^{1/3}. The unit of a_{c} can be fixed by multiplying the r.h.s. of this proportionality with R^{2/3}. Otherwise the proportionality coefficient can only depend on μ_{T}, and of course, on the sign of Δγ. Therefore, we can write
such that the r.h.s. of the equation only depends on the Tabor coefficient and the functional form of the surface interaction. Since we have not used the explicit functional form of our default surface interaction (other than that it depends only on a single length scale), the conclusions drawn in this section extends to any choice for V_{fr} considered in this work.
To include finite loads into the analysis, note that the ratio F_{N}/ΔγR does not change under the transformation Equation 13–Equation 17. This allows us to express a properly undimensionalized contact radius as a dimensionless function of a properly dimensionless load
From this last equation, it also becomes clear that the pulloff (or the squeezeout) force is proportional to ΔγR, i.e., by identifying the value of F_{N}/ΔγR at which the function a_{c,T} takes its minimum value. Therefore, it is most meaningful to normalize the force with ΔγR, unless, of course, Δγ = 0.
The approach is validated in Figure 3. It shows the spatial dependence of the gap for two different parameterizations. Small deviations, which are not visible to the naked eye, occur. They are due to finitesize and discretization effects. For example, the ratio a_{c}/L is not exactly zero and takes different values for different values of s for a fixed number of points used to represent the elastic surface.
The normal displacement can be nondimensionalized in a similar fashion as the contact radius, except that it needs to be rescaled with the factor w^{2/3} rather than with w^{1/3}. This is why it must obey
where all terms on the r.h.s. of the equation are again dimensionless, while those on the l.h.s. are allowed to have units. Thus, displacements and gaps are best represented in units of R^{1/3}(Δγ/E)^{2/3}, while contact radii are more meaningfully expressed as multiples of R^{2/3}(Δγ/E)^{1/3}. As a result, numbers turn out of order unity when data is represented in these units, unless F_{N} approaches the pulloff threshold or distinctly exceeds ΔγR.
I conclude this section by summarizing the units used in this work and discuss some of the consequences arising from it: Specifically, the following units are used for:
This list includes a “new” unit of normal stress or pressure, [p], which must be chosen proportional to [F_{N}]/[x]^{2} rather than to E so that the regular definition of pressure applies. As noted above, E drops out of the definition of the unit for the normal force, implying that pulloff or squeezeout forces cannot be functions of E. Instead they must equal RΔγ times dimensionless expressions that can only depend on μ_{T}, and, of course, on the functional form of the interaction potential. In our units, the wellknown a_{c}(F_{N},E^{*},R,Δγ) relations can be simplified to
Hertzian contact mechanics is obtained in either limit for F_{N} >> 1. Finally, note that Maugis’ choice for units slightly differs from ours in that he used πΔγ rather than Δγ in Equation 23–Equation 26 and 3E/4 instead of E. The conversion between Maugis’ and our system is summarized in Equation 43–Equation 46.
Numerical analysis
Different methods can be used in the numerical solution of Equation 11. For the present study, Green’s function molecular dynamics (GFMD) [24] as described in [25] was employed. The only modification is the implementation of conservative surface forces acting in addition to the boundary condition g(r) ≥ 0. Moreover, the results in this work were produced with a serial code with typical run times of a few minutes. I refer to the literature for more details on GFMD [24,25]. Irrespective of the employed code or method, particular precautions, which are worth discussing, need to be taken into account when including adhesion or finiterange repulsion.
When simulating Hertzian contact mechanics, one needs to ensure that the discretization of the lattice Δa satisfies Δa << a_{c}. Of course, methods based on spatially varying grids only need to obey that relation near the contact line. In addition, one wants a_{c} to be much less than the size of the simulation cell, at least in Fourierbased methods, such as GFMD. One then has the sequence of inequalities Δa << a_{c} << L. In Hertzian contact mechanics, this is easy to achieve: choosing the discretization such that Δa/a_{c} = 1/32 and a_{c}/L = 1/8 already gives accurate results for the contact area, that is, to within less than 0.1% error if the contact area is determined through a fit of the radial pressure profile.
When including adhesion, an additional length enters, namely that associated with the adhesive zone. The additional adhesive radius or skin a_{a} then needs to be taken into consideration. When the Tabor coefficient is very small, a_{a} becomes large, and one needs a_{a} to lie within the simulation cell. A new series of inequalities is obtained: Δa << a_{c} << a_{a} << L. If the normal stress changes smoothly with the gap, i.e., for longrange adhesion, the forces couple predominantly to large wavelength modes. This then results in a simple offset force, as is well known from DMT theory. However, numerical demands can become significant when a_{c} disappears continuously with decreasing load as in the DMT scenario. The condition a_{c} << a_{a} then becomes difficult to satisfy.
In the opposite case of a large Tabor coefficient, a_{a} is very small, potentially much smaller than a_{c}. We still need to resolve the adhesive zone sufficiently well, because the stress has to be smooth on the given discretization. One thus obtains the series of inequalities a << a_{a} << a_{c} << L. In either limit of large or small Tabor coefficient, another inequality needs to be satisfied in addition to those for Hertzian contact mechanics.
While the contact area converges reasonably quickly as the respective inequalities are obeyed, the centerofmass displacement d, which corresponds to or u(r → ∞) only converges slowly. The reason is that the displacement field induced at a given point due to an external force only decays with the inverse distance from that point, i.e.,
where c is a load and systemdependent constant. Outside of the adhesive zone, this relation can be used, in principle, to extrapolate from finite r to infinite r, i.e., by determining c and d from two measurements taken sufficiently far outside the adhesive zone. In practice, this turns out problematic, because the periodic boundary conditions suppress the 1/r corrections near the boundaries.
For Hertzian contacts, it is still possible to use a slightly modified (empirical) correction
where X and M denote the (symmetry) points (L/2,L/2) and (L/2,0) relative to the position of the center of the tip. The same extrapolation scheme also appears to give quickly converging estimates for the normal displacement for adhesive contact, which is demonstrated in Figure 4.
It is worth discussing Figure 4. At the given load, the contact radius is a_{c} ≈ 2.30, while it would have been a_{c} ≈ 1.05 without adhesion. The displacement curve has a peak at r = 2.52 and adhesive effects remain nonnegligible all the way up to r ≈ 4. At that distance the gap starts to be greater than 5z_{0}, which means that the adhesive attraction is less than e^{−5} times the value in the contact and its immediate periphery. For distances exceeding r = 4, an infinitely large system would then show the displacement given in Equation 29. The periodic boundary conditions suppress this scaling rather strongly, yet, for radii as small as r = 10, accurate estimates for d_{∞} can be achieved through Equation 30.
Simulations could be made more efficiently by exploiting the radial symmetry of the system. This would allow one to reduce sums over two indices (e.g., q_{1} and q_{2}) to that over one index. However, less is gained than it first might seem. To get a good resolution of the contact area, the onedimensional (1D) calculation require greater ratios for a_{c}/a than twodimensional setups. The reason is that the resolution of the contact area in 1D and in 2D both scale with 1/N, where N is the number of grid points in the contact. For example, when representing a contact in which for the given discretization 5a < a_{c} < 6a in a 2D system, then a_{c} is allowed to take the values and The maximum distance between two radii thus is Δa_{max} = 0.28 so that the resolution is Δa_{c}/a ≈ 5/0.28. To match this in a 1D model, one would need 18 grid points rather than five. Since we need to cover a (square) area of (2a_{c})^{2} in 2D, we thus have a computational overhead of a little more than a factor of 4 compared to 1D. However, the disadvantage of large 1D systems is that the number of iterations to find solutions can be much larger than in 2D. Depending on the nature of the solver, the number of iterations scales as a power law with linear system size. In the given case, where the effective stiffness scales proportional to wave vector q, one would expect a slowing down with at least in simple gradientbased minimization such as steepest descent or moleculardynamics. For this reason, no efforts were made to reduce the dimensionality, although this would have been legitimate for the given problem.
Positive work of adhesion
This section analyzes how the employed models reproduce established results for adhesive singleasperity contacts in the limits of large and small Tabor coefficient. This includes an asymptotic analysis of Maugis’ solution of the Dugdale model, which in turn leads to modifications of the equations proposed by Carpick, Ogletree, and Salmeron [1]. The crossover from JKR to DMT is investigated as well, in particular at zero load and near pulloff, allowing one to identify the model for the surface interaction that is most appropriate for the simulation of (adhesive) multiasperity contacts.
Zero external load
An external load of F_{N} = 0 is addressed first. The motivation for studying this special case is that one can analyze relatively easily at what Tabor coefficients the DMT and JKR limits start to predict reasonably accurate values for the contact radii and displacements in our various models.
We start our analysis with the pressure distribution of the exponential model, which is depicted in Figure 5 for μ_{T} = 1/4 and μ_{T} = 4. It behaves very similar to the MD model, which shall not be shown explicitly. As to be expected, the adhesive load is spread out for μ_{T} = 1/4 to radii clearly exceeding a_{c} (all the more as each radius r contributes with a weight proportional to r), while it is rather localized near r = a_{c} for μ_{T} = 4. It therefore seems legitimate to call the (net) pressure profile for μ_{T} = 1/4 DMTlike and JKRlike for μ_{T} = 4.
The adhesive pressure is calculated from the functional derivative p_{adh}(x,y) = −δV_{fr}/δg(x,y), where V_{fr} is defined in Equation 5. This can be evaluated to yield
which becomes p_{adh}(r < a_{c}) = −Δγ/z_{0} within the true contact area. Using Equation 19, one obtains
Stress or pressure originating from the constraint g(x,y) ≥ 0 is computed from the elastic Green’s functions.
The wellknown qualitative difference for the contact geometry of systems with large (μ_{T} = 4) and small (μ_{T} = 1/4) Tabor coefficient is borne out in the radial dependence of the gap g(r). Specifically, Figure 6 reveals that a small Tabor coefficient makes g(r) have a positive curvature at r ≥ a_{c}, as in the DMT solution, while it has a negative curvature at r ≥ a_{c}, indicative of an adhesive neck, for large μ_{T}. Figure 6 also shows that the displacement (defined by the difference between the actual gap and the gap in an undeformed contact at r >> 1) is smaller for μ_{T} = 4 than for μ_{T} = 1/4, although the contact radius is larger for μ_{T} = 4 than for μ_{T} = 1/4.
The Gauss model behaves qualitatively different from the exponential model. First, there are no adhesive forces within the contact, but only outside of it, as shown in Figure 7. Second, at r = a_{c}, the total pressure disappears in the Gauss model, while it remains finite in the exponential model. Third, the pressure due to the constraint has finite slope at r = in the Gauss model, while the slope diverges in the exponential model (not shown explicitly). All these differences arise because the derivative of v_{fr}(g/z_{0}) remains finite (i.e., with positive sign) in the exponential model when the gap closes, while is zero in the Gauss model.
Another consequence of p_{adh}(g/z_{0}) having zero slope in the limit of g → 0^{+} is that the gap in the Gauss model closes continuously rather than with a discontinuity in its first derivative. This is shown in Figure 8, from where it becomes clear that it is difficult to define good measures for the contact radius. In a linear representation (at low resolution), it seems as if the contact closes with the typical adhesive neck, i.e., in part (a) of Figure 8 the gap appears to close at r ≈ 2.415. There, the slope of g(r) takes its maximum value, in a very similar fashion as in the JKR limit, or for the exponential model for the same value of μ_{T} = 16. However, when increasing the magnification, one can see that the contact closes only at r ≈ 1.97. Unfortunately, the radii where the gap closes to zero, and where g(r) has its maximum do not approach each other quickly when μ_{T} is increased. Instead, the value of g in the crossover region in Figure 8b moves to smaller values as μ_{T} increases. Similar behavior is seen in the VDW model, which is not shown here explicitly.
Zeroload contact radii for different potentials are depicted in Figure 9 as function of the Tabor coefficient. In the exponential model, the contact radius approaches DMT and JKR limits in a very similar fashion as in the MD model. In a later section on the asymptotic analysis, I find that the MD corrections to the JKR limit are of order for large Tabor coefficients while those to the DMT limit are of order μ_{T} for small Tabor coefficients. The same scaling of the leadingorder corrections is apparently borne out in the exponential model.
Models in which v_{fr}(g/z_{0}) has zero slope in the origin behave qualitatively different from the MD or the exponential model. They approach the DMT limit for the contact radius fairly quickly, i.e., roughly with However, convergence to the JKR limit is poor. The latter can be improved by defining the contact line to be located where g′(r) takes its maximum value. Unfortunately, this definition cannot be universally applied, i.e., only when μ_{T} is sufficiently large to allow for an adhesive neck to be formed, see also Figure 8. Moreover, in the context of randomly rough surfaces with complicated contact geometries, this last definition of contact would not be practicable.
Unlike the contact radius, the normal displacement d does not suffer from any difficulties to be properly defined. In principle, this could enable one to ascertain v_{fr}(g) from displacement measurements without much ambiguity. However, Figure 10 reveals that the functional form of d(μ_{T}) is relatively insensitive to the details of the finiterange interaction, at least, as long as we allow for a redefinition of the Tabor coefficient, such that all curves superimpose at the distance half way between the JKR and the DMT limit. This is in agreement with a work by Tvergaard and Hutchinson [18] who found that Δγ and the peak stress (which one may losely associate with Δγ/z_{0}) are the basic parameters for mode I fracture.
Before proceeding to the case of finite load, I wish to comment on the relatively large numerical (GFMD) uncertainties for the displacement at large Tabor coefficients. They stem predominantly from the difficulty to apply the finitesize extrapolation formula, Equation 30, to gaps having adhesive necks. This problem would not be present in largescale simulations of multiasperity interfaces, because system sizes would automatically be much larger than local contact radii. One may conclude that the use of the exponential model for the study of adhesive multiasperity contacts appears to be appropriate. The MD model could be used as well, in principle, however, it might induce undesired numerical instabilities due to the discontinuity of v′(g/z_{0}) at g = z_{0}. The Gauss model can only be taken when the property of interest is related to the gap but not for the calculation of contact area. If one wanted to simulate van der Waals attraction at large μ_{T}, one might want to replace the VDW model in Equation 12 with 1/(1 + g/z_{0})^{2}. This dependence makes it possible to determine the contact area meaningfully in the realm of continuum mechanics while using reasonable approximations for van der Waals interactions at large distance.
Finite external load
In most experiments, the Tabor coefficient is kept constant and the normal load is changed. As a result, one obtains the normal displacement d(F_{N}) as a function of the normal load F_{N}. In some cases, i.e., for sufficiently large contact radii, an estimate of the contact radius, a_{c}(F_{N}), can be obtained as well. One might be tempted to believe that knowing such curves allows one to deduce the surface forces. Here, I want to investigate to what degree such an inversion is possible by studying the sensitivity of the functions d(F_{N}) and a_{c}(F_{N}) to the functional form of the surface interactions. Figure 11 shows the contact radius as a function of the normal load. One can see that the exponential model and the MD model agree very closely, that is, curves almost superimpose for a given Tabor coefficient. This makes it essentially impossible to discriminate between these two forms of interaction experimentally. Likewise, the Gauss and VDW models also coincide for the same Tabor coefficient despite their significant differences at large gaps. Interestingly, the μ_{T} = 1 curve for VDW and Gauss (both having finite slope potentials at g = 0) is akin of the μ_{T} = 1/4 curves for the MD and the exponential model (both having zeroslope potentials at g = 0). This confirms the trend reflected in Figure 9: Surface potentials with zero slope at g = 0 make the results move toward the DMT limit.
In contrast to the a_{c}(F_{N}) dependence, the normal displacement curve d(F_{N}) predominantly depends on the Tabor coefficient. Now, all μ_{T} = 1 curves resemble each other closely, independent of the slope of the surface potential at zero gap. As for the normal displacement, all curves are reasonably close to the JKR limit. Even the μ_{T} = 1/4 curve lies closer to the JKR than to the DMT line. This is consistent with the results shown in Figure 10, which show that the DMT limit for d(F_{N}) is only reached at extremely small Tabor coefficients. Figure 12 reveals that it is possible to adjust the free parameters of the MD model to fit d(F_{N}) curves for a broad variety of surface interactions. However, one should abstain from deducing contact areas based on such fits, as this can result in nonnegligible errors. For example, if we only knew the contact area from Maugis’ solution, we would be illadvised to conclude from Figure 12 that the contact area for the μ_{T} = 1 Gauss model should lie roughly half way between those of the μ_{T} = 1 and μ_{T} = 4.
Comparison to other models and asymptotic analysis
Maugis proposed an analytical solution for the relation between contact radius a_{c} and normal load F_{N} in the Dugdale model. It requires the elimination of an auxiliary variable, m, through the selfconsistent solution of two coupled nonlinear equations. Once a_{c} and m are found, the displacement d can be readily calculated as well. Using tildes to indicate variables in Maugis’ unit system, the relevant equations read:
where the functions f(m), g(m), h(m), and j(m) are defined as
and
In each but one (straightforward) case, behavior of the functions for m approaching unity or infinity has been included. They become useful in the limit of large and small Tabor coefficients, respectively.
Conversion back to our unit system can be done using:
To overcome the need of having to find the selfconsistent solution to Maugis’ equations, Carpick, Ogletree, and Salmeron (COS) [1] proposed a simple and thus elegant analytical formula for the a_{c}(F_{N}) dependence
Schwarz [2] later recognized that the COS description is exact – given proper parameterization – if the interaction between the surfaces results from the superposition of an infinitesimally shortranged and an infinitely longranged contribution. However, in the given context of one intermediaterange potential, I will treat the COS equation as a guessed approximation containing the correct functional form in the limits of large and small μ_{T}.
The primary COS equation (Equation 47) is designed such that the contact radius at zero load a_{c}(0,μ_{T}) as well as the pulloff force F_{p}(μ_{T}) can be reproduced exactly. However, approximations to their dependence on μ_{T} had been provided as well, because no closedform expression are available. A free parameter remains, α(μ_{T}), which can be used to minimize deviations from the exact solution. At large loads, one recovers the wellknown a_{c} scaling, however, not necessarily with the correct prefactor. Another property of the COS approximation is that it does not necessarily contain the correct value of the contact radius at pulloff. Thus, despite predicting the contact radius fairly well, the COS contact radius is not exact in the limit of very large and very small (i.e., pulloff) normal loads. These deficiencies can be improved when parameterizing the COS equation in a slightly different fashion, e.g.,
with
This set of equation ensures that a_{c} converges to the exact value when F_{N} → ∞ and F_{N} → −F_{p}. The parameter α(μ_{T}) can then be adjusted to either yield the correct zeroload contact radius, or to minimize the deviation between approximation and the exact Maugis solution by some other mean. Note that the factors 3/4 in Equation 48 and 4/3 in Equation 49 have to be replaced by unity when working with Maugis’ unit system.
I wish to note that including the correct asymptotics in the a_{c}(F_{N}) expression does not necessarily improve the fits in the range from slightly above the pulloff force at negative loads to several times the absolute pulloff force. This is demonstrated in Figure 13. Moreover, convergence to the correct a_{c}(F_{N}) dependence at large loads is rather slow even when using Equation 49.
It is also possible to constrain the COS relation for the contact radius such that it contains the correct pulloff force and contact radius as well the correct zeroload radius. In either case, relative errors are small, i.e., ≤1% even for μ_{T} ≈ 1, where one is relatively far away from both the DMT and the JKR limit.
Zero load: The asymptotic analysis is readily done for zero loads, because the variable m can be directly eliminated in that case. As a result, one obtains
and
From the last two equations, one can see – as in Figure 9 and Figure 10 – that the JKR limit is quickly reached as the Tabor coefficient increases. However, convergence to the DMT limit with decreasing is rather slow. It is particularly slow for the normal displacement. E.g., to have a maximum error in a_{c}(F_{N} = 0) and d(F_{N} = 0) that is of order 1% with respect to a desired limit, it is sufficient to work with ≈ 10 for JKR, but one needs ≈ 10^{−4} to approach the DMT limit with similar accuracy. The latter is not problematic for the simulation of multiasperity contacts, as the system is large by default. However, for singleasperity contacts, large deviations from μ_{T} = 1 (on a logarithmic scale) are difficult to handle in singleasperity contact simulations for reasons discussed in the numericalanalysis section.
Knowing the asymptotic behavior of and with respect to allows one to incorporate it in empirical equations for these two quantities. The following equations are found to achieve this and to provide excellent approximations to Maugis’ solutions:
Two coefficients in each of the last two equations (c_{1}, c_{2} and c_{5}, c_{6}) can be constrained to reproduce the correct asymptotics (and thus be obtained analytically). Two fit parameters remain for contact radius and one for the displacement. The relative errors from the pertinent fits are shown in Figure 14. Compared to an already quite accurate empirical relation proposed by Carpick et al. for , see Eq. (12b) in [1], the new Equation 52 and Equation 53 contain the correct asymptotics and reduces the maximum relative error from 1.5% to 0.3%. The data shown in Figure 14 were obtained with the following numerical values: c_{1} = 4/5, c_{2} = −1.285, c_{3} = 4/5, c_{4} = −0.435, c_{5} = −3/2, c_{6} = 0.1845, and c_{7} = 6.71.
Asymptotic behavior near pull off: The structure of the COS approximation, Equation 47, and its modified form, Equation 48, indicates that the critical behavior near pull off satisfies where κ must be either κ = 1/3 as in the DMT limit, or κ = 1/2 as in the JKR limit. However, nothing in the selfconsistent solution of Maugis indicates that the exponent κ changes discontinuously from one value to the next as the Tabor coefficient reaches or passes through a critical value. In fact, representing the data from Figure 11 in terms of Δa_{c}(ΔF_{N}), as done in Figure 15, shows that κ changes continuously from 1/3 to 1/2 as μ_{T} increases from 0 to infinity.
An analysis for the normal displacement, similar to the one presented in Figure 15 but not shown explicitly here, exhibits a similar trend. The exponent describing Δd = d − d_{p} as a function of ΔF_{N} = F_{N} + F_{p} crosses over continuously from the DMT to the JKR limit as μ_{T} increases. However, there is not a onetoone relation between μ_{T} and κ. In particular the data sets for μ_{T} = 1 show relatively large differences between the exponent in the MD model (κ ≈ 0.469) and the exponential model (κ ≈ 0.435).
The insights obtained from Figure 15 can be used, in principle, to further modify the COS approximations, e.g., by replacing the squareroot in Equation 48 by some other power or likewise by changing the squareroot and the exponent 2/3 on the r.h.s. of Equation 47 in an appropriate fashion. When doing so, the modified version of Equation 48 does not only converge to the correct value at pulloff. It can also be parameterized to yield the correct asymptotics near pulloff. This results in a further reduction of the mean or overall error by a little more than a factor of two with respect to those shown in Figure 13, however, at the expense of one additional fit variable. Since the main new aspect of this study is concerned with negative work of adhesion and, moreover, both original and modified COS equations are already quite accurate, a more detailed analysis of the adhesive singleasperity contact is not pursued in this work.
Negative work of adhesion
For repulsive contacts, Δγ < 0, there is obviously no finite contact radius at zero normal load F_{N} = 0^{+}. The repelled rigid tip simply “hovers” at (infinitely) large distance over an undeformed elastic manifold. This is why it is not possible in this case to conduct a zeroload analysis similar to that presented for adhesive contacts. Since Maugis’ solution has not yet been extended to repulsive contacts, we are not in a position to compare our data to analytical solutions for negative Δγ. One of the consequences is that the asymptotic analysis must be based on GFMD data, except for μ_{T} → 0, for which normal forces couple predominantly to longwavelength modes so that the Hertzplusoffset approximation (DMT) should be accurate. Given the close similarity between the exponential and the Maugis–Dugdale model as well as that between the Gauss and the VDW model seen in the last section, the attention is restricted to one potential in each class, i.e., the exponential and the Gauss model.
We start our analysis with the contact radius dependence on load. In analogy to the context of wetting fluids, one may call the force at which a finite value of a_{c} becomes unstable upon lowering the load the spontaneous wetting force F_{sw}. The force above which a_{c} can no longer be zero is called the squeezeout force F_{sq}. If the transition from contact to noncontact is continuous F_{sw} = F_{sq}, otherwise F_{sw} < F_{sq}. Results are shown in Figure 16.
As is the case for attractive interactions, the contact radius at small loads can be sensitive to both the Tabor coefficient and the choice of the potential. Specifically, the exponential model always shows a continuous transition from finite to zero contact radius (at least for the values of μ_{T} investigated here), while the Gauss model has either a continuous transition below a critical Tabor coefficient ≤ 1 or a discontinuous transition for μ_{T} > The discontinuity of the contact radius for Gauss potentials and sufficiently large Tabor coefficients implies that two solutions may coexist, i.e., one where the two surfaces are separated and one where they touch. However, once F_{N} exceeds a second threshold force F_{sq}(μ_{T}), i.e., the squeezeout force, only one solution survives, that is, the one with finite contact radius. This can be seen in analogy to adhesive contacts with μ_{T} > 0, where two solutions coexist in a finite interval of forces −F_{p} ≤ F_{N} ≤ 0.
As for the a_{c}(F_{N}) relation near pulloff in the case of positive work of adhesion, the excess contact radius, a_{c} − a_{sw}, depends as a power law on the excess force, F_{N} − F_{sw}, for F_{N} ≥ F_{sw}:
Fits to the a_{c}(F_{N}) relation are shown exemplarily for two values of μ_{T} in Figure 17. Details about the fits to the presented as well as additional data are summarized in Table 1. As for attractive contacts, it is found that κ changes continuously from κ(μ_{T} → ∞) = 1/2 to κ(μ_{T} → 0) = 1/3. For small μ_{T} , Hertzplusoffset behavior is reached as evidenced by the observation that c and F_{sw} approach (3/4)^{1/3} and 2π, respectively. However, F_{sw} as well as F_{sq} quickly increase with μ_{T} for μ_{T} ≥ 1. This latter behavior is different from that of the pulloff F_{p} force for attractive surfaces, which only varies between 1.5π and 2π in the present unit system. Since the increase of both F_{sw} and F_{sq} is much faster in the exponential model than in the Gauss model, one can conclude that the exponential model converges more quickly to the continuum model than the Gauss model.
Table 1: Results of fits to the data shown in Figure 16. The last digit may not be significant.
model  μ_{T}  a_{sw}  F_{sw}  F_{sq}  c  κ 

Gauss  1/16  0  6.30  6.30  0.91  0.333 
1/4  0.01  6.43  6.44  0.86  0.34  
1  0.25  8.00  8.40  0.56  0.46  
4  0.66  47.3  86  0.30  0.50  
exponential  1/16  0  6.38  6.38  0.62  0.35 
1/4  0  6.89  6.89  0.68  0.44  
1  0  13.85  13.85  0.46  0.48  
4  0  339  339  0.28  0.49 
As in the case of adhesive interactions, the normal displacement seems less sensitive to both the choice of the potential and the Tabor coefficient than the contact area, unless normal loads are very small, i.e., at loads similar in magnitude or smaller than the squeezeout load for μ_{T} = 1. This is demonstrated in Figure 18. It reveals that information on the (effective) nearrange surface interactions at small separation are difficult to obtain from experimentally measured loaddisplacement curves.
I conclude this section with an analysis of the gap profile for repulsive contacts. At large loads, different Tabor parameters and functional forms for finiterange repulsion yield gap profiles that are indistinguishable at small magnification, see Figure 19a. Differences become nevertheless significant at high resolution near the center of the contact. Particularly remarkable is the data set for the Gauss model with μ_{T} = 4 and its bistability revealed in Figure 19b. For an increasing force, no contact has formed at F_{N} = 7.5. However, when reaching F_{N} = 7.5 from above, contact is formed for radii r < a_{c} ≈ 1.73. In the latter case, the gap then quickly increases within Δr ≈ 0.1 to an almost constant value of order 1/μ_{T} for r ≥ a_{c}, as if one had a single confined layer of liquid. For radii r > , the gap assumes the “macroscopic” behavior. Here, ≈ 4 is the contact radius that one would ascertain from the analysis of the gap with low resolution, e.g., via graphical inspection of Figure 19a.
At small loads, the sensitivity of the gap profile on the details of the model become even more apparent. This result, which can be seen in Figure 19c, is expected, since the elasticity of the tip is no longer relevant. Instead, the forcedisplacement curve is predominantly determined by the effective surface interactions, as shown clearly by the μ_{T} = 4 data sets in Figure 18. They exhibit, to leading order, a F_{N} exp(−d/ζ) relation in the exponential model, and a F_{N} exp(−d^{2}/ζ^{2}) relation for the Gauss model, where ζ is inversely proportional to μ_{T}.
Conclusion
The principle new aspect of this work is the continuummechanics based analysis of singleasperity contacts with finiterange repulsion acting in addition to shortrange hardwall repulsion. The analysis is based on the concept of the Tabor coefficient and the repulsion is assumed to arise due to the presence of a strongly wetting fluid. As for attractive singleasperity contacts, it is found that the contact area or the displacement on the normal load depend, to a large degree, not only on the surface energy but also on the Tabor coefficient μ_{T}. Moreover, for μ_{T} exceeding a critical value, there may exist a range of loads in which two (meta)stable solutions coexist, i.e., one in which the surfaces touch and one in which a thin gap between the two surfaces remains. When the value for the load is increased above a threshold, the latter solution becomes unstable and the gap disappears. However, in order to obtain this kind of behavior, which is reminiscent of the squeezeout of a wetting fluid, the finiterange interactions between the contacting surfaces have to be tailored correctly. Using a surface interaction v_{fr}, whose derivative increases monotonically as the gap g approaches zero, such as v_{fr} exp(−g/z_{0}), only one stable solution exists for any given normal load. Conversely, when the distance–force dependence is multivalued, as is the case for a v_{fr} relation, squeezeout and spontaneous wetting can be rationalized and thus be modeled in the realm of continuum mechanics – in terms of transitions between (meta)stable solutions. These transitions (similar to instabilities in the Prandtl model [26], in which a particle is dragged with a weak spring through a sinusoidal potential) can occur for solvated tips on surfaces, for example, if the effective tip–surface interactions has zero slope when the surfaces touch, as is the case for v_{fr} . In reality, the farfield potential may even be oscillatory and evidenced by the squeezeout of many subsequent layers. Such behavior has been recently observed and linked to the (damped) longrange oscillatory behavior of the density correlations in highdensity liquids [15,19].
An interesting consequence of shortrange repulsion is that the contact geometry can look similar to that of an adhesive neck. This is shown in Figure 19b for the (μ_{T} = 4) Gauss model and decreasing load. To improve the visualization, a similar gap geometry is shown again in Figure 20 together with a profile of the finiterange repulsion.
A secondary aspect of this work is devoted to the analysis of how to best reach welldefined asymptotic behavior in numerical simulations of adhesive contact mechanics. It is found that the DMT limit is approached quickest when using attractive potentials whose first derivative disappears as the gap goes to zero, at least if the contact area is the variable of interest. However, these potentials approach the JKR limit only at a rate of 1/μ_{T} for large μ_{T} and the contact area becomes difficult to define once μ_{T} ≥ 1. Thus, one is better off using potentials with finite slope in the smallgap limit. They converge in a welldefined fashion with to the JKR limit for large Tabor coefficients. This is supposedly the more relevant limit for adhesive surfaces with selfaffine fractal roughness. For the modeling of repulsive surfaces, the situation is more complicated. Formally, the JKR limit is again reached more quickly with models that have finite slope at zero gap. However, these models do not allow one to model the hysteretic response of a confined fluid that results whenever the squeezeout force exceeds the spontaneous wetting force.
A byproduct of this work is a minor modification of the phenomenological description of singleasperity contact mechanics by Carpick, Ogletree, and Salmeron [1]. The COS equations can be parametrized to contain the correct asymptotic behavior for JKR and for DMT limits and also for the superposition of extremely short and longrange interfacial interactions, as shown by Schwarz [2]. However, they still have a few formal shortcomings for intermediaterange potentials. For example, the original interpolation of the contactareaonload dependence for finite Tabor coefficients recuperates neither Hertzian contact mechanics at large loads with correct prefactors nor the correct contact radii in either DMT or JKR limit at zero normal loads. In this work, I propose to enforce those limits exactly including the correct asymptotics for a_{c}(µ_{T},F_{N} = 0) at μ_{T} = 0 and μ_{T} → ∞. By doing so, the maximum error of the a_{c}(F_{N} = 0,μ_{T}) curve could be reduced from 1.2% to less than 0.3%. A shortcoming of both the original and the new, modified COS equations is that they both assume an asymptotic behavior near pulloff (F_{N} → −F_{p}) according to (a − a_{p}) (F_{p} + F_{N})^{κ}, where the exponent takes the JKR value κ = 2/3 for any nonzero Tabor coefficient. The modified COS equations could thus be improved further if one incorporated the new finding that κ crosses over continuously from 2/3 (exact for JKR) to 1/2 (exact for DMT). However, this does not seem useful in practice. Extreme accuracies (5 digits and more for a_{c} and F_{N}) would be needed in measurements to deduce a_{p} and κ to within one or two digits. Such an accuracy is difficult to achieve both experimentally and numerically. Moreover, the surface energy is not very well defined at small scales, because its precise value depends crucially on roughness down to the atomic scale, see, e.g., [27]. Thus, from a practical point of view, both the original and the modified COS equations are quite reasonable, all the more because the geometry of real tips can deviate quite substantially from a parabola.
This work is concluded with an assessment of what values for μ_{T} one might expect in AFM or SFA experiments. To come up with a ballpark estimate, the following “typical values” shall be assumed: Δγ = 40 mN (Δγ can, of course, be close to zero, but much higher, e.g., for two equally charged surfaces in the context of electrochemistry), E = 5 GPa (in between soft matter and ceramics), z_{0} = 10 Å (size of an OMCTS or molten salt molecule), R = 1 μm (in between AFM and SFA, precise value not very important, as third root is taken). These numbers lead to μ_{T} = 0.4, which is close to the interesting “crossover” regime. Thus, real contacts may span a broad range of values for μ_{T}. Comparison between theory and experiment may be difficult, in particular because atomicscale roughness (or even subatomic roughness arising from electron orbitals) leads to complicated slipboundary conditions and slow kinetics. However, given a wellmotivated form for the effective interaction between two flat surfaces, it may yet be possible to rationalize and to model, at least on a semiquantitative level, the interactions of curved surfaces in the presence of a strongly wetting fluid within the presented Taborcoefficient based framework. Particularly appealing systems may be found in triboelectrochemical applications, where the surface interactions can be tailored in a quasicontinuous fashion.
References

Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027
Return to citation in text: [1] [2] [3] [4] [5] [6] [7] 
Schwarz, U. D. J. Colloid Interface Sci. 2003, 261, 99–106. doi:10.1016/S00219797(03)000493
Return to citation in text: [1] [2] [3] 
Grierson, D. S.; Flater, E. E.; Carpick, R. W. J. Adhes. Sci. Technol. 2005, 19, 291–311. doi:10.1163/1568561054352685
Return to citation in text: [1] 
Hertz, G. J. Reine Angew. Math. 1881, 92, 156. doi:10.1515/crll.1882.92.156
Return to citation in text: [1] 
Derjaguin, B. V.; Muller, V. M.; Toporov, Yu. P. J. Colloid Interface Sci. 1975, 53, 314–326. doi:10.1016/00219797(75)900181
Return to citation in text: [1] 
Johnson, K. L.; Kendall, K.; Roberts, A. D. Proc. R. Soc. London, Ser. A 1971, 324, 301–313. doi:10.1098/rspa.1971.0141
Return to citation in text: [1] 
Tabor, D. J. Colloid Interface Sci. 1977, 58, 2–13. doi:10.1016/00219797(77)903666
Return to citation in text: [1] 
Muller, V. M.; Yushenko, V. S.; Derjaguin, B. V. J. Colloid Interface Sci. 1980, 77, 91–101. doi:10.1016/00219797(80)904191
Return to citation in text: [1] 
Maugis, D. J. Colloid Interface Sci. 1992, 150, 243–269. doi:10.1016/00219797(92)90285T
Return to citation in text: [1] [2] 
Barthel, E. J. Phys. D: Appl. Phys. 2008, 41, 163001. doi:10.1088/00223727/41/16/163001
Return to citation in text: [1] 
Hughes, B. D.; White, L. R. Q. J. Mech. Appl. Math. 1979, 32, 445–471. doi:10.1093/qjmam/32.4.445
Return to citation in text: [1] 
Vinogradova, O. I.; Feuillebois, F. J. Colloid Interface Sci. 2003, 268, 464–475. doi:10.1016/j.jcis.2003.09.002
Return to citation in text: [1] [2] 
Persson, B. N. J.; Tosatti, E. J. Chem. Phys. 2001, 115, 5597–5610. doi:10.1063/1.1398300
Return to citation in text: [1] 
Müser, M. H. Phys. Rev. Lett. 2008, 100, 055504. doi:10.1103/PhysRevLett.100.055504
Return to citation in text: [1] 
Fisher, M. E.; Wiodm, B. J. Chem. Phys. 1969, 50, 3756–3772. doi:10.1063/1.1671624
Return to citation in text: [1] [2] [3] [4] 
Chandra, N.; Li, H.; Shet, C.; Ghonem, H. Int. J. Solids Struct. 2002, 39, 2827–2855. doi:10.1016/S00207683(02)00149X
Return to citation in text: [1] 
Turon, A.; Dáviall, C.; Camanho, P.; Costa, J. Eng. Fract. Mech. 2007, 74, 1665–1682. doi:10.1016/j.engfracmech.2006.08.025
Return to citation in text: [1] 
Tvergaard, V.; Hutchinson, J. W. Int. J. Solids Struct. 1996, 33, 3297–3308. doi:10.1016/00207683(95)002618
Return to citation in text: [1] [2] 
Hoth, J.; Hausen, F.; Müser, M. H.; Bennewitz, R. J. Phys.: Condens. Matter 2014.
accepted.
Return to citation in text: [1] [2] 
Luan, B. Q.; Robbins, M. O. Nature 2005, 435, 929–932. doi:10.1038/nature03700
Return to citation in text: [1] 
Mo, Y. F.; Turner, K. T.; Szlufarska, I. Nature 2009, 457, 1116–1119. doi:10.1038/nature07748
Return to citation in text: [1] 
Shengfeng, C.; Robbins, M. O. Tribol. Lett. 2010, 39, 329–348. doi:10.1007/s1124901096825
Return to citation in text: [1] 
Eder, S.; Vernes, A.; Vorlaufer, G.; Betz, G. J. Phys.: Condens. Matter 2011, 23, 175004. doi:10.1088/09538984/23/17/175004
Return to citation in text: [1] 
Campañá, C.; Müser, M. H. Phys. Rev. B 2006, 74, 075420. doi:10.1103/PhysRevB.74.075420
Return to citation in text: [1] [2] 
Dapp, W. B.; Lücke, A.; Persson, B. N. J.; Müser, M. H. Phys. Rev. Lett. 2012, 108, 244301. doi:10.1103/PhysRevLett.108.244301
Return to citation in text: [1] [2] 
Prandtl, L. Z. Angew. Math. Mech. 1928, 8, 85–106. doi:10.1002/zamm.19280080202
Return to citation in text: [1] 
Jacobs, T. D. B.; Ryan, K. E.; Keating, P. L.; Grierson, D. S.; Lefever, J. A.; Turner, K. T.; Harrison, J. A.; Carpick, R. W. Tribol. Lett. 2013, 50, 81–93. doi:10.1007/s1124901200973
Return to citation in text: [1]
18.  Tvergaard, V.; Hutchinson, J. W. Int. J. Solids Struct. 1996, 33, 3297–3308. doi:10.1016/00207683(95)002618 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
2.  Schwarz, U. D. J. Colloid Interface Sci. 2003, 261, 99–106. doi:10.1016/S00219797(03)000493 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
2.  Schwarz, U. D. J. Colloid Interface Sci. 2003, 261, 99–106. doi:10.1016/S00219797(03)000493 
3.  Grierson, D. S.; Flater, E. E.; Carpick, R. W. J. Adhes. Sci. Technol. 2005, 19, 291–311. doi:10.1163/1568561054352685 
7.  Tabor, D. J. Colloid Interface Sci. 1977, 58, 2–13. doi:10.1016/00219797(77)903666 
16.  Chandra, N.; Li, H.; Shet, C.; Ghonem, H. Int. J. Solids Struct. 2002, 39, 2827–2855. doi:10.1016/S00207683(02)00149X 
6.  Johnson, K. L.; Kendall, K.; Roberts, A. D. Proc. R. Soc. London, Ser. A 1971, 324, 301–313. doi:10.1098/rspa.1971.0141 
17.  Turon, A.; Dáviall, C.; Camanho, P.; Costa, J. Eng. Fract. Mech. 2007, 74, 1665–1682. doi:10.1016/j.engfracmech.2006.08.025 
5.  Derjaguin, B. V.; Muller, V. M.; Toporov, Yu. P. J. Colloid Interface Sci. 1975, 53, 314–326. doi:10.1016/00219797(75)900181 
14.  Müser, M. H. Phys. Rev. Lett. 2008, 100, 055504. doi:10.1103/PhysRevLett.100.055504 
2.  Schwarz, U. D. J. Colloid Interface Sci. 2003, 261, 99–106. doi:10.1016/S00219797(03)000493 
15.  Fisher, M. E.; Wiodm, B. J. Chem. Phys. 1969, 50, 3756–3772. doi:10.1063/1.1671624 
27.  Jacobs, T. D. B.; Ryan, K. E.; Keating, P. L.; Grierson, D. S.; Lefever, J. A.; Turner, K. T.; Harrison, J. A.; Carpick, R. W. Tribol. Lett. 2013, 50, 81–93. doi:10.1007/s1124901200973 
11.  Hughes, B. D.; White, L. R. Q. J. Mech. Appl. Math. 1979, 32, 445–471. doi:10.1093/qjmam/32.4.445 
12.  Vinogradova, O. I.; Feuillebois, F. J. Colloid Interface Sci. 2003, 268, 464–475. doi:10.1016/j.jcis.2003.09.002 
13.  Persson, B. N. J.; Tosatti, E. J. Chem. Phys. 2001, 115, 5597–5610. doi:10.1063/1.1398300 
15.  Fisher, M. E.; Wiodm, B. J. Chem. Phys. 1969, 50, 3756–3772. doi:10.1063/1.1671624 
19. 
Hoth, J.; Hausen, F.; Müser, M. H.; Bennewitz, R. J. Phys.: Condens. Matter 2014.
accepted. 
10.  Barthel, E. J. Phys. D: Appl. Phys. 2008, 41, 163001. doi:10.1088/00223727/41/16/163001 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
9.  Maugis, D. J. Colloid Interface Sci. 1992, 150, 243–269. doi:10.1016/00219797(92)90285T 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
8.  Muller, V. M.; Yushenko, V. S.; Derjaguin, B. V. J. Colloid Interface Sci. 1980, 77, 91–101. doi:10.1016/00219797(80)904191 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
26.  Prandtl, L. Z. Angew. Math. Mech. 1928, 8, 85–106. doi:10.1002/zamm.19280080202 
15.  Fisher, M. E.; Wiodm, B. J. Chem. Phys. 1969, 50, 3756–3772. doi:10.1063/1.1671624 
18.  Tvergaard, V.; Hutchinson, J. W. Int. J. Solids Struct. 1996, 33, 3297–3308. doi:10.1016/00207683(95)002618 
15.  Fisher, M. E.; Wiodm, B. J. Chem. Phys. 1969, 50, 3756–3772. doi:10.1063/1.1671624 
24.  Campañá, C.; Müser, M. H. Phys. Rev. B 2006, 74, 075420. doi:10.1103/PhysRevB.74.075420 
25.  Dapp, W. B.; Lücke, A.; Persson, B. N. J.; Müser, M. H. Phys. Rev. Lett. 2012, 108, 244301. doi:10.1103/PhysRevLett.108.244301 
1.  Carpick, R. W.; Ogletree, D. F.; Salmeron, M. J. Colloid Interface Sci. 1999, 211, 395–400. doi:10.1006/jcis.1998.6027 
24.  Campañá, C.; Müser, M. H. Phys. Rev. B 2006, 74, 075420. doi:10.1103/PhysRevB.74.075420 
25.  Dapp, W. B.; Lücke, A.; Persson, B. N. J.; Müser, M. H. Phys. Rev. Lett. 2012, 108, 244301. doi:10.1103/PhysRevLett.108.244301 
12.  Vinogradova, O. I.; Feuillebois, F. J. Colloid Interface Sci. 2003, 268, 464–475. doi:10.1016/j.jcis.2003.09.002 
9.  Maugis, D. J. Colloid Interface Sci. 1992, 150, 243–269. doi:10.1016/00219797(92)90285T 
19. 
Hoth, J.; Hausen, F.; Müser, M. H.; Bennewitz, R. J. Phys.: Condens. Matter 2014.
accepted. 
20.  Luan, B. Q.; Robbins, M. O. Nature 2005, 435, 929–932. doi:10.1038/nature03700 
21.  Mo, Y. F.; Turner, K. T.; Szlufarska, I. Nature 2009, 457, 1116–1119. doi:10.1038/nature07748 
22.  Shengfeng, C.; Robbins, M. O. Tribol. Lett. 2010, 39, 329–348. doi:10.1007/s1124901096825 
23.  Eder, S.; Vernes, A.; Vorlaufer, G.; Betz, G. J. Phys.: Condens. Matter 2011, 23, 175004. doi:10.1088/09538984/23/17/175004 
© 2014 Müser; licensee BeilsteinInstitut.
This is an Open Access article under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
The license is subject to the Beilstein Journal of Nanotechnology terms and conditions: (http://www.beilsteinjournals.org/bjnano)