Distribution of Gaps and Adhesive Interaction Between Contacting Rough Surfaces

Understanding the distribution of interfacial separations between contacting rough surfaces is integral for providing quantitative estimates for adhesive forces between them. Assuming non-adhesive, frictionless contact of self-affine surfaces, we derive the distribution of separations between surfaces near the contact edge. The distribution exhibits a power-law divergence for small gaps, and we use numerical simulations with fine resolution to confirm the scaling. The characteristic length scale over which the power-law regime persists is given by the product of the rms surface slope and the mean diameter of contacting regions. We show that these results remain valid for weakly adhesive contacts and connect these observations to recent theories for adhesion between rough surfaces.


Introduction
Contact between nominally flat, rough surfaces has been the subject of study for countless experimental, analytical, and numerical investigations over the past century . A universal theme found in most cases is that contact is limited to the peaks or asperities of the rough topography and the real area of contact A rep is much less than the apparent projected area A 0 . Substantial progress has been made in determining the relationship between A rep and the applied normal force F in non-adhesive, frictionless systems assuming elastic response. In such cases, the proportionality A rep = rep F∕h � 0 E * is found when A rep ≲ 0.1A 0 . Here, h ′ 0 is the root mean square (rms) slope of the rough topography, E * is the elastic contact modulus, and the dimensionless constant rep ≈ 2 [2,6,7]. This result implies a load-independent mean compressive stress, rep = h � 0 E * ∕ rep in contacting regions. Since flattening a region with slope h ′ 0 introduces a strain of order h ′ 0 into the surface, rep reflects the stress required to flatten the rough topography.
Recently, interest has turned to systems including attractive interactions that lead to macroscopic adhesion. Two opposite limits exist in the classical literature on the adhesion of smooth spheres: in the Derjaguin-Muller-Toporov (DMT) [32] limit, weak attractive forces between solids do not alter the geometry of contact, but do reduce the global mean pressure; in the opposite limit, known as the Johnson-Kendall-Roberts (JKR) [33] limit, strong attractive forces significantly change the structure of the contact edge and can lead to contact hysteresis. For the contact of spheres, the Tabor and Maugis parameters [34] describe the continuous transition between the two limits in terms of the relative strengths of adhesive and elastic parameters.
Working in a DMT-like limit, Pastewka and Robbins developed a theory to predict the onset of stickiness in contacts of self-affine rough surfaces [23] that was later independently confirmed by Müser [35]. They split the total normal force as F = F rep + F att into the sum of a repulsive contribution F rep > 0 and an attractive contribution F att < 0 . Repulsive and attractive contributions originate from repulsive surface patches of total area A rep and attractive surface patches of area A att (see Fig. 1). The total force is then given by where att is the mean stress in the attractive patches that is roughly constant and of order att = w∕Δr , where w is the work of adhesion and Δr the range of the attractive interaction.
Like in the purely non-adhesive limit, the geometry of contact is fractal in the DMT-like limit with proportionality between A rep and the contact perimeter P rep given by where the mean contact diameter d rep is approximately constant [7,23,36]. Note that Eq. (2) holds generally for any geometric object, but d rep varies with A rep in most cases. Short-ranged attractive interactions generate narrow bands of approximately constant width d att located around contact regions (see Fig. 1). If d att is small, then the total area contributing to attractive forces A att = P rep d att , which can be related to A rep via mutual proportionality with P rep given by Eq. (2). This means, that for non-sticky interfaces, the (repulsive) contact area is given by the expression with an effective 1∕ = 1∕ rep − 1∕ att . The adhesive interaction hence increases the effective value of the dimensionless constant . A macroscopic force is required to separate the two surface when |F att | ≈ F rep or equivalently att ≈ rep . Interfaces that require a macroscopic force for separation are called "sticky".
This theory depends sensitively upon the distribution of interfacial separations or gaps, that is assumed to be unaltered from the non-adhesive scenario. Previous work has primarily focused on the behavior of the mean gap ḡ , which is commonly found to be exponentially related to the normal load in non-adhesive contact as F ∝ exp −ḡ∕ h 0 , where h 0 is the rms surface height and is a dimensionless constant of order unity [8, 9, 12-17, 21, 24, 29]. Almqvist et al. [16], the contact mechanics challenge [29] and Wang and Müser [37] have reported distributions of interfacial separations, but these works have either not focused on the behavior at small gaps that is important for understanding short-ranged adhesion or indirectly reported it through analysis of percolation in Reynolds flow.
In this paper, we derive the distribution of interfacial separations in the vicinity of contacting regions and show numerically that our expression holds even in the weakly adhesive limit. This distribution can be used to compute the total attractive contribution to the force and hence the force-area relationship for weakly adhesive interfaces.

Simulation Methods
In our simulations, we invoke the standard mapping that allows the contact of two rough, elastic solids to be treated as contact between an initially flat, elastic solid and a rough, rigid surface [38]. If the elastic properties of the two original surfaces are encoded by the Young's moduli E 1 and E 2 and Poisson's ratios The roughness profile of each periodic, L × L surface with nominal area A 0 = L 2 is described by a self-affine fractal between an upper cutoff length scale max and a lower cutoff min , where length scales are given in terms of the pixel size a 0 . This means that the power spectral density (PSD) C( ) of the isotropic, self-affine roughness depends only on wavevector magnitude q = | | and satisfies where q 1 = 2 ∕ max and q 2 = 2 ∕ min are the wavevector magnitudes corresponding to the roughness cutoffs and H is the Hurst exponent ( 0 < H < 1 ) that determines correlations in the roughness. The resolution of the calculations is given by min ∕a 0 , while the ratio L∕ max controls the surface "representativity" [18]. The distribution of heights P(h) for a Gaussian self-affine rough surface with mean height h and L∕ max ≫ 1 has (by construction) Gaussian form. We use min ≥ 32a 0 for the results presented here to ensure that the contact edge is sufficiently resolved. We use a Fourierfiltering algorithm to create such self-affine surfaces for our calculations [36,39].
Performing the mapping described above results in a rough surface that is the incoherent sum of the rough profiles of the original surfaces. For the combination of two profiles that have identical statistical properties, the commonly utilized statistical measures of roughness-the rms height We therefore here work in the limit of a rigid rough surface contacting a flat deformable elastic halfspace and all statistical properties are to be interpreted for a combined surface. Note that ⟨⋅⟩ denotes the spatial average over the domain where the topography function h(x, y) is defined.
We use a static boundary element method to compute the linear elastic deformation induced by normal contact on an isotropic half-space (Refs. [40][41][42]). In order for linear elasticity to be a good approximation, the surface slope must be small. We use h � 0 = 0.1 to ensure that this approximation is justified, but note that real surfaces may have values h � 0 ≈ 1 or larger [31,43,44].
Our simulations assume frictionless contact with a noninterpenetration constraint. Pixels on the surface of the elastic solid are considered to be in contact when they bear a compressive pressure; each contacting pixel contributes an area a 2 0 to the total contact area. For non-adhesive calculations, we solve the interpenetration constraint using a constrained conjugate-gradient optimizer [45]. For adhesive calculations, we assume an interaction energy of that depends on the gap g, where is the interaction range. (Note that the range defined in Ref. [23] is Δr ≈ 1.36 for this potential.) The overall adhesive energy is then given by where g(x, y) is the local gap at position x, y. This attractive interaction is minimized with an interpenetration constraint realized through the constrained non-linear conjugate-gradient algorithm of Ref. [46]. Similar boundary element methods have been used extensively to study rough contacts [11, 18-23, 25-28, 30]. The attractive interaction is identical to the one employed in the contact mechanics challenge [29].

Theory
The distribution of interfacial separations g between a rough surface and an undeformed elastic solid with surface at h = 0 is equal to the (Gaussian) distribution of heights, where ḡ is the initial mean surface height. The width of the distribution is the same as for the rough surface itself and scales with the long wavelength cutoff, h 0 ∼ H max . When the solids are pushed together under load, the elastic solid deforms and the mean interfacial separation ḡ shrinks. We now write the gap distribution as the additive decomposition where p c (g) contains the distribution of the gaps g within the contacting area, p n (g) the contribution from near the contact edge, and p f (g) the contribution from farther distances from the contact edge. Since p(g) is normalized, the probability of contact for a given contact fraction c = A rep ∕A 0 is where is the Dirac -function.
The next contribution p n (g) comes from small separations near the contact edge. As shown in numerical calculations in Refs. [7,23], the contact edge has a total perimeter P rep = A rep ∕d rep with a constant d rep . (This means that both area and contact edge are fractal objects-see Fig. 1a-and suggests that their fractal dimensions are identical.) We can write this contribution as where Δ(x) describes how the mean gap between the two contacting surfaces varies as a function of distance x from the contact edge (see Fig. 1b). Note that we take the integral in Eq. (11) out to a characteristic length that is close g(x, y)), p(g) = p c (g) + p n (g) + p f (g), enough to the contact edge such that the number of points contributing to the gap distribution is still proportional to P rep . The arguments that lead to Eq. (11) rely on the geometry of the contact patches formed between contacting self-affine surfaces. Refs. [7,23] showed that in this case P rep ∝ A rep , compared to P rep ∝ A 1∕2 rep as expected for simple contact shapes like circles. The larger perimeter scaling exponent arises because fractal contact regions are not compact, and a significant perimeter contribution comes from regions inside the convex hull enclosing individual patches. As patches become larger, they simply contain more non-contacting regions. The presence of non-contacting regions makes the contact geometry appear locally rectangular. The deformation cross-section for each line segment drawn through a single continuous contact diameter looks like that of a cylindrical indenter, rather than a spherical one. Within Euclidean geometry, a rectangle with constant thickness along its minor axis has the property P rep ∝ A rep , if the rectangle is thin enough such that the contributions to P rep from the shorter sides can be neglected. The characteristic width of the contacting rectangle, and hence the contact diameter for the cylinder, is d rep .
Working with this analogy, we now derive an analytic prediction for the distribution of gaps produced by nonadhesive contact of smooth surfaces. Assuming a non-adhesive cylindrical contact [47,48], we have where h ′ 0 is the local slope of the cylinder at the contact edge. It can be shown generally for non-adhesive contact that the local separation Δ(x) ∝ x 3∕2 for small lateral distances x from the contact edge [38]. The characteristic scale for the interfacial separation at the contact edge is g 0 = 4h � 0 d rep ∕3 , the prefactor in Eq. (12). Inserting Eqs. (12) into (11) yields where Eq. (2) was used for the right-hand equality. This is our prediction for the gap distribution near the contact edge. For small g, the distribution diverges as g −1∕3 .
We note that ∫ ∞ 0 dg p(g) = 1 , but taking the integral out to infinity for the contribution to p(g) given by Eq. (13) diverges. The contribution p n (g) to the gap distribution can therefore only be valid up to g ∼ g 0 , the characteristic gap that is reached a distance d rep from the contact edge. The "far" contribution p f (g) looks like the undeformed distribution of heights given by Eq. (8) (see also Refs. [16,29]).
Almqvist et al. [16] have observed a divergence of the gap distribution for small g but have not quantified the exponent.
Pastewka and Robbins [23] have used this divergence for their theory of "stickiness" but have not provided extensive numerical evidence for its validity. We now supplement these observations with additional high-resolution numerical data.

Results
We performed simulations with sufficiently fine resolution ( min = 32a 0 and 128a 0 ) to test Eq. (13). Figure 2 shows distributions of interfacial separations for non-contacting grid points. The distributions were normalized to the respective non-contacting fractional area, 1 − c , and divided by the prefactor of Eq. (13) to collapse all data points onto a single (g∕g 0 ) −1∕3 power-law. We used the numerically measured value of d rep (see Ref. [23] on details of how to compute d rep ) for the data collapse. The power-law regime emerges for both H = 0.3 (Fig. 2a) and H = 0.8 (Fig. 2b). The data only collapse over a limited range of gaps. The divergence at small gap predicted by Eq. (13) is always cut off at a minimum length scale, below which p(g) is uniformly distributed. In our simulations, the cutoff scale is ∼ 10 −3 − 10 −2 a 0 ; the threshold for saturation of p(g) at small g∕g 0 is inversely proportional to the resolution a 0 ∕ min . For H = 0.3 (Fig. 2a), contributions to p(g) from near-contact and far-from-contact overlap (i.e., p n (g) and p f (g) have similar magnitude) for min = 32a 0 (solid symbols). The power-law regime only clearly emerges over 1-2 decades for H = 0.3 if the resolution of the calculation is improved by increasing min , as is shown for min = 128a 0 (open symbols). For H = 0.8 , on the other hand, the power-law regime extends over a much larger range of gaps even for our "coarse" calculations with min = 32a 0 . As for H = 0.3 , increasing the short-wavelength cutoff to min = 128a 0 extends the power-law to smaller gaps.
One can integrate p n (g) to calculate the cumulative fractional area c n (g c ) in the near-contacting region closer than a cutoff gap g c . We find which for g c = g 0 yields c n ∕c = . This means the area in the near-contact region is ∼ 3 times the area in the contacting regions. The area in the near-contact region is equal to the contact area when g c ≈ 0.18g 0 . The breakdown of the power-law region at small g occurs at g∕g 0 ≈ 10 −3 or smaller (for H = 0.8 ), meaning that the area within the rolloff region at small g is at most about 3% of the contact area.
For large gaps the power-law regime is cut off by a Gaussian gap distribution that reflects the distribution of undeformed or weakly deformed parts of the surface (see also (14) Eq. (8)). For both H = 0.3 and H = 0.8 , the uptick in p(g) in the range g∕g 0 ∼ 1 − 10 is the peak of this Gaussian height distribution. Since h 0 ∝ ( max ∕ min ) H , the power-law regime extends much further for H = 0.8 (Fig. 2b) than for H = 0.3 (Fig. 2a). Increasing max ∕ min shifts the Gaussian peak out to larger g∕g 0 (most prominently for H > 0.5 ), and extends the range over which p(g) ≈ p n (g) . However, for the largest max ∕ min we studied, p(g) < p n (g) for 0.1 ≲ g∕g 0 ≲ 1 . This suggests that the cutoff of the integral in Eq. (11) may be up to an order of magnitude smaller than g 0 , about equal to the point where the near-contact area matches the contact area. Still larger calculations may be required to verify this result.
The crucial question for adhesive theories is how applicable our results are for gap distributions in the presence of attractive interactions. We quantify the strength of the attractive interaction by the value of 1∕ att (see Ref. [23] and discussion below), since interfaces become sticky for 1∕ att ≳ 1∕2 [23]. For the data collapse, we use the values of d rep measured in the non-adhesive calculations. Figure 3 shows the gap distributions for adhesive calculations as w increases (with constant ), using surfaces with H = 0.8 and min = 128a 0 . To facilitate the comparison between non-adhesive and adhesive calculations, the adhesive simulations are conducted with the constraint that the mean gap is identical to the non-adhesive simulation ( 1∕ att = 0 ). This choice of constraint means that the contact areas are not equal, particularly for sticky surfaces. For non-sticky surfaces (up to 1∕ att ≈ 0.2 ), the gap distribution follows the predicted power-law over the same range as the non-adhesive result. At 1∕ att ≈ 0.2 ( A rep ∕A 0 ≈ 0.06 ) there is a slight deviation toward smaller gaps. The sticky cases with 1∕ att ≈ 1 and 1∕ att ≈ 2 ( A rep ∕A 0 ≈ 0.12 and 0.16, respectively) clearly deviate from our prediction for the divergence. These sticky interfaces appear to still exhibit a regime where p(g) ∝ g −1∕3 , but the range of the power-law regime becomes narrower as 1∕ att increases. This means that the prefactor from Eq. (13) no longer captures the intensity of the divergence and much of the non-contacting area is pushed out toward larger gaps. This is also reflected by the increase of the peak at larger gaps. The characteristic gap below which the distribution rolls-off and becomes constant appears to increase slightly in the sticky limit.

Discussion and Conclusions
The behavior of the near-contact interfacial separation is an important consideration in the context of adhesive contact because it determines the attractive contribution to the force, Assuming that the contribution from p f (g) is negligible, i.e., that our potential is sufficiently short-ranged, the attractive pressure per unit contact area is then given by where we have used the interaction law v(g) from Eq. (6). Note that Eq. (16) defines the value of the dimensionless constant 1∕ att = F att ∕h � 0 E * A rep , which we used to quantify the strength of adhesion in Fig. 3 and which can be used to determine the effective range of adhesion Δr (see Ref. [23]).
This expression can also be used to quantify what "short-ranged" adhesion means for rough surfaces and, therefore, to determine the limits of the theories of Pastewka and Robbins [23] and Müser [35]. Figure 4 shows the normalized integrand of Eq. (16) as a function of the normalized gap g∕g 0 . The integrand depends on the range of the interaction potential . For ≈ 0.1g 0 and below, the main contribution to the integrand and thereby the attractive force comes from gaps with g∕g 0 < 0.1 where the power-law holds. The roll-off region at small gap contributes negligibly to this integral. This means that the adhesion range is short enough if ≲ 0.1g 0 . The calculations presented here were carried out with = 4a 0 (as is typical for attractive interactions, e.g., Refs. [49,50]) and our adhesive calculations have g 0 > 25a 0 , which means the range is sufficiently small. Interactions with a larger , range (e.g., electrostatic interactions) interact with the full topography of the surface and require corrections to the expression for F att . The same is true for extremely smooth surfaces with small values of g 0 . This defines bounds for theory outlined in Refs. [23,35]. Besides leaving the gap distribution unmodified, the DMT-like limit also implies that the repulsive area A rep equals the typical definition of the contact area, namely, vanishing gaps g = 0 . The latter definition makes sense in continuum theories (like the present work), but not for models that consider the full intermolecular interaction, which exhibits soft repulsion (e.g., Ref. [23]), or for those that include thermal fluctuations (e.g., Refs. [51,52]). Like in the JKR model [33], these definitions of contact no longer agree for sticky interfaces. In addition, JKR-like contacts separate as x 1∕2 and not x 3∕2 near the contact edge [53], leading to a gap distribution p n (g) ∝ g that is clearly distinct from what we observe in our calculations. Furthermore, in the sticky limit, the contact geometry changes substantially and is not expected to give rise to the simple proportional contact area-perimeter relationship used in our calculations. This is especially true for soft solids, which can deform to fill in interior non-contacting regions without large elastic energy penalties, thereby increasing the contact area at the expense of the overall contact perimeter.
Nevertheless, the present calculations are a powerful demonstration of the universal emergence of the g −1∕3 divergence in the distribution of gaps between elastically stiff rough surfaces. Since this behavior is a direct consequence of the fractal character of the contacting interfaces as manifested in the proportionality between perimeter and contact area, our results are another indirect demonstration of this aspect of the contact geometry. The distribution is unaltered by weak adhesive interactions, giving additional support for Fig. 4 Integrand of Eq. (16), i.e., non-dimensionalized contribution to the overall attractive force as a function of gap g∕g 0 for different ranges of the adhesive interaction the DMT-like approximation that underlies the adhesive theories of Pastewka and Robbins [23] and Müser [35].