Single-asperity contact mechanics with positive and negative work of adhesion: Influence of finite-range interactions and a continuum description for the squeeze-out of wetting fluids

Summary In this work, single-asperity contact mechanics is investigated for positive and negative work of adhesion Δγ. In the latter case, finite-range repulsion acts in addition to hard-wall 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 load-bearing or a squeezed-out fluid. The possibility for coexistence and the subsequent discontinuous wetting and squeeze-out instabilities depend not only on the Tabor coefficient μT but also on the functional form of the finite-range repulsion. For example, coexistence and discontinuous wetting or squeeze-out 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 finite-range attraction. The results can benefit the interpretation of atomic force microscopy in liquid environments and the modeling of multi-asperity contacts.


Introduction
The continuum description of single-asperity contact mechanics has received much attention in the last few decades. This is in large parts because it describes force-displacement curves rather accurately down to nanometer scales relevant to atomic force microscopy (AFM) [1][2][3]. The contributions to the linear elasticity of (frictionless) single-asperity contacts most central to this work are the following: Hertz [4] solved the contact mechanics of a parabolic tip pressed against a substrate for hard-wall 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 long-range interaction, in which case adhesion effectively acts as a normal load that is independent of the contact-radius a c . Johnson, Kendall, and Roberts (JKR) [6] solved the problem for short-range 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 longand short-range 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 cohesive-zone model introduced by Dugdale (MD) and found analytical solutions for intermediate-range adhesion at arbitrary values of μ T .
Although single-asperity, linearly-elastic, 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 finite-range repulsion between two surfaces acting in addition to hard-wall 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 squeeze-out as well as a comparison to the behavior near pull-off 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 pull-off force and pull-off 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 squeeze-out (finite-range repulsion) versus pull-off (finite-range 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 contact-mechanics between macroscopic, adhesive, rough surfaces in the context of continuum-mechanics, 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 finite-range 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 single-asperity 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 single-asperity contact problem is introduced. As shown in Figure 1, we consider a stiff, ideally-flat wall positioned in the xy plane at z = 0 and a linearly-elastic tip of parabolic shape. Its undeformed surface is given by (1) where R is the radius of curvature and the in-plane 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 in-plane 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 non-holonomic boundary condition (3) Alternatively, one can formally introduce a short-range, hardwall repulsion [14] (4) 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 energy-based description of the short-range repulsion. Figure 1: Geometry of the deformed tip (upper grey solid), the substrate (lower solid), and the reference tip (solid line). The dotted line represents a hypothetical tip that is allowed to penetrate the substrate the distance d into the substrate without deforming. The following vectors are introduced in the sketch: Normal load F N acting on the center of mass of the tip, the elastic displacement field u, and the displacement d of the tip's center of mass. In addition, two scalar quantities, namely the contact radius a c and the gap (field) g are shown.
This work also considers finite-range adhesive or finite-range repulsive interactions V fr , which only depend on the local gap. The default expression for it is: (5) 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 intermediate-range interactions can be modeled as well. Note that V fr and V sr are qualitatively different: The prefactor of the short-range 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 finite-range 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 in-plane 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 non-negligible adhesive (or finite-range 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 small-slope approximation) as (8) 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 (9) When solving the contact problem, one seeks to minimize the total energy (10) with respect to u(x,y), i.e., the solution u 0 (r) must satisfy (11) 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 traction-separation 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: (12) 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 long-range density correlations in fluids are governed by a single length [15]. In a high-density fluid, the correlation length becomes complex [15], which then leads to layering transitions as discussed recently [19]. The VDW model might approximate the long-range 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 hard-wall repulsion) and that attraction -or additional repulsion -may occur at finite distance. Continuous short-range 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 hard-wall or hard-disk interactions [20][21][22][23]. Lastly, the equations to be solved would become very stiff and thus the simulation inefficient if the hard-wall constraint was replaced by potentials with large curvature.
In the context of the squeeze-out 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 squeeze-out 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 finite-range 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 in-plane 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 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 hard-wall repulsive energy V fr is unproblematic. It disappears for the old and the new solution, because of the limit z r → 0, i.e., .
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 long-ranged can only depend on a non-dimensionalized 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 (18) 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 (19) 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 in-plane 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 (20) 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 (21) From this last equation, it also becomes clear that the pull-off (or the squeeze-out) 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 finite-size 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 (22) 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 pull-off 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 pull-off or squeeze-out 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 well-known a c (F N ,E * ,R,Δγ) relations can be simplified to (27) (28) 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

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 finite-range 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 Fourier-based 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 long-range 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 center-of-mass 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 system-dependent 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 (30) 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.  Figure 1) of a tip with R = E = 1 pressed against a rigid substrate at a positive load of F N = 1.5625 for a work of adhesion Δγ = 1 and a Tabor coefficient μ T = 1 resulting in an exponential decay length of z 0 = 1. In each case, the system is discretized into 512 × 512 grid points, but different sizes are used, i.e., L = L 0 , L = 2 1/3 L 0 , L = 2 2/3 L 0 with L 0 = 9.8825. Part (a) shows a larger domain including the shape of the tip in form of a dashed line. Part (b) shows a smaller domain and includes a higherresolution estimate of the displacement at infinite radius R through the extrapolation 6u X − 5u M .
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 non-negligible 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 one-dimensional (1D) calculation require greater ratios for a c /a than two-dimensional 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 gradient-based minimization such as steepest descent or molecular-dynamics. 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 single-asperity 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 cross-over from JKR to DMT is investigated as well, in particular at zero load and near pull-off, allowing one to identify the model for the surface interaction that is most appropriate for the simulation of (adhesive) multi-asperity 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 DMT-like and JKR-like 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 (31)  Stress or pressure originating from the constraint g(x,y) ≥ 0 is computed from the elastic Green's functions.
The well-known 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 cross-over 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. Zero-load 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 leading-order 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 finite-range 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 finite-size extrapolation formula, Equation 30, to gaps having adhesive necks. This problem would not be present in large-scale simulations of multi-asperity 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 multi-asperity 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 exponen- Figure 11: Contact radius a c as a function of load F N for the exponential and the MD model using different Tabor coefficients, ranging from μ T → ∞ (JKR, top) to μ T = 0 (DMT, bottom). For the Gauss and VDW models, only μ T = 1 is shown. Their a c (F N ) curve is similar to that of the Maugis and the exponential model for μ T = 1/4. Color coding: μ T = 4 (red), μ T = 1 (green), and μ T = 1/4 (blue). tial 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 zero-slope 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 non-negligible errors. For example, if we only knew the contact area from Maugis' solution, we would be ill-advised 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.  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.

Comparison to other models and asymptotic analysis
Conversion back to our unit system can be done using: To overcome the need of having to find the self-consistent 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 (47) Schwarz [2] later recognized that the COS description is exactgiven proper parameterization -if the interaction between the surfaces results from the superposition of an infinitesimally short-ranged and an infinitely long-ranged contribution. However, in the given context of one intermediate-range 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 pull-off force F p (μ T ) can be reproduced exactly. However, approximations to their dependence on μ T had been provided as well, because no closed-form 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 well-known 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 pull-off. 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., pull-off) normal loads. These deficiencies can be improved when parameterizing the COS equation in a slightly different fashion, e.g., 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 pull-off force at negative loads to several times the absolute pull-off 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 pull-off force and contact radius as well the correct zero-load radius. In either case, rela- tive 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 (50) and (51) 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 multi-asperity contacts, as the system is large by default. However, for single-asperity contacts, large deviations from μ T = 1 (on a logarithmic scale) are difficult to handle in single-asperity contact simulations for reasons discussed in the numerical-analysis 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.  where κ must be either κ = 1/3 as in the DMT limit, or κ = 1/2 as in the JKR limit. However, nothing in the self-consistent 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. Here a c and a p denote the contact radius at an arbitrary load F N and the pull-off load F p , respectively. Deviations between the Maugis solution and the exponential model are particularly obvious for the μ T = 1 data set. Color coding: μ T = 4 (red), μ T = 1 (green), and μ T = 1/4 (blue).
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 one-toone 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 square-root in Equation 48 by some other power or likewise by changing the square-root 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 pull-off. It can also be parameterized to yield the correct asymptotics near pull-off. 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 single-asperity 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 zero-load 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 long-wavelength modes so that the Hertz-plus-offset 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 squeeze-out force F sq . If the transition from contact to non-contact 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 squeeze-out force, only one solution survives, that is, the one with finite contact radius. This can be seen in analogy to adhesive contacts with   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 , Hertz-plusoffset 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 pull-off 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. 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 squeeze-out 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 load-displacement 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 finite-range 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 Figure 18: Displacement d as a function of normal force F N in the vicinity of the spontaneous wetting force F sw . Symbols reflect numerical results. The lines, which connect many data points not explicitly shown, are drawn to guide the eye. The two thick grey lines reflect the square of the contact radius in the Hertz and DMT approximation, respectively. Color coding: μ T = 4 (red), μ T = 1 (green), and μ T = 1/4 (blue).

Figure 19:
Gap g(r) as a function of the lateral distance from the origin r for a large load F N = 60 (a) and (b) as well as for an intermediate load F N = 7.5 (c). In each case, surfaces repel each other. Graph (b) contains the same data as (a) but has higher resolution. Color coding: μ T = 4 (red) and μ T = 1 (green). 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 force-displacement 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 single-asperity contacts with finiterange repulsion acting in addition to short-range hard-wall 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 single-asperity 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 squeeze-out of a wetting fluid, the finite-range 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, squeeze-out 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 far-field potential may even be oscillatory and evidenced by the squeeze-out of many subsequent layers. Such behavior has been recently observed and linked to the (damped) long-range oscillatory behavior of the density correlations in high-density liquids [15,19].
An interesting consequence of short-range 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 finite-range repulsion.
A secondary aspect of this work is devoted to the analysis of how to best reach well-defined 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 small-gap 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 self-affine 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 squeeze-out force exceeds the spontaneous wetting force.
A by-product of this work is a minor modification of the phenomenological description of single-asperity 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 superpo-sition of extremely short-and long-range interfacial interactions, as shown by Schwarz [2]. However, they still have a few formal shortcomings for intermediate-range potentials. For example, the original interpolation of the contact-area-on-load 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 pull-off (F N → −F p ) according to (a − a p ) (F p + F N ) κ , where the exponent takes the JKR value κ = 2/3 for any non-zero 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 "cross-over" regime. Thus, real contacts may span a broad range of values for μ T . Comparison between theory and experiment may be difficult, in particular because atomic-scale roughness (or even sub-atomic roughness arising from electron orbitals) leads to complicated slip-boundary conditions and slow kinetics. However, given a well-motivated form for the effective interaction between two flat surfaces, it may yet be possible to rationalize and to model, at least on a semi-quantitative level, the interactions of curved surfaces in the presence of a strongly wetting fluid within the presented Tabor-coefficient based framework. Particularly appealing systems may be found in tribo-electrochemical applications, where the surface interactions can be tailored in a quasi-continuous fashion.