---
author:
- 'Vladislav A. Yastrebov'
title: |
    The Elastic Contact of Rough Spheres\
    Investigated Using a Deterministic Multi-Asperity Model
---

In this paper we use a deterministic multi-asperity model to investigate the elastic contact of rough spheres. Synthetic rough surfaces with controllable spectra were used to identify individual asperities, their locations and curvatures. The deterministic analysis enables to capture both particular deformation modes of individual rough surfaces and also statistical deformation regimes, which involve averaging over a big number of roughness realizations. Two regimes of contact area growth were identified: the Hertzian regime at light loads at the scale of a single asperity, and the linear regime at higher loads involving multiple contacting asperities. The transition between the regimes occurs at the load which depends on the second and the fourth spectral moments. It is shown that at light indentation the radius of circumference delimiting the contact area is always considerably larger than Hertzian contact radius. Therefore, it suggests that there is no scale separation in contact problems at light loads. In particular, the geometrical shape cannot be considered separately from the surface roughness at least for approaching greater than one standard roughness deviation.

[**Keywords.**]{} roughness; indentation; contact; deterministic multi-asperity model; contact area.

Introduction
============

Contact and friction interactions play an essential role in many quotidian contexts, including those related to industry and transportation (e.g., tire-road and wheel-rail contacts, electric switches, gears, bearings, and brake systems), everyday human activity (e.g., walking, handling, touching, and sitting) and natural phenomena (e.g., earthquakes, landslides, and glacier motion). Regardless of such prevalence, contact-related mechanisms (friction, adhesion, and wear) are still not fully understood and thus are among the most cutting-edge research topics in the mechanical community [@mo2009friction; @ben2010dynamics; @pastewka2011anisotropic; @aghababaei2016critical].

Numerous models of contact-related mechanisms exist at structural scale [@vakis2018]. They serve to model (1) interfacial normal and tangential stiffness [@wriggers1993application; @akarapu2011stiffness], (2) frictional resistance [@ruina1983slip], (3) material removal on rubbing surfaces (wear) [@archard1956wear; @yin2010adhesive], (4) heat transfer between contacting solids [@madhusudana1996thermal], (5) contact electric resistance [@slade2017electrical; @yastrebov2015holm], (6) adhesion [@fuller1975effect; @pastewka2014contact], (7) interfacial fluid flow [@dapp2012self; @dapp2016fluid; @shvarts2018trapped], (8) microstructural changes in near-contact material layers [@ramesh2008modeling], (9) fretting wear life-cycle [@hills1994mechanics; @fouvry1997wear; @dick2006fretting], (10) debris generation [@sayles1982influence; @basseville2011numerical], (11) lubrication [@patir1978average], and other mechanisms. The associated models can be incorporated in a macroscopic/structural model via constitutive interfacial equations. These equations can be based either on experimental data, and thus remain purely phenomenological, or can take the microscopic roughness as the starting point. The latter class of models shall have a greater predictive power, and potentially can be used for a large spectrum of applications. However, because of the strong non-linearity of the contact/friction mechanisms and extreme complexity of surface roughness, construction of a reliable analytical micro-mechanical model presents a serious challenge. Moreover, the scale separation which would allow to separate the macroscopic shape of contacting solids and the accompanying roughness cannot be always ensured.

Admittedly, all the aforementioned phenomena are strongly related to the surface roughness. (1) The compression of contacting asperities as well as their interlocking contributes to the interfacial stiffness. (2) According to the adhesive theory of friction the frictional resistance of contacting spots is pressure independent and thus the total frictional resistance of the interface is proportional to the true contact area. In general, this area may depend on the hold-time or sliding velocity manifesting itself as a velocity-dependent friction coefficient. (3) In case of adhesive wear, the removed volume depends on the characteristics and number of contact spots forming the true contact area. (4) The overall thermal conductivity of the interface is determined by heat (phonons and, if available, electrons) passing directly through contact spots, by convective heat transfer through the interface fluid and also by the radiative heat transfer. All these heat exchange mechanisms depend on the true contact area and on the opening (or gap) field between non-contacting parts of the two solids. (5) In a rather similar way, the electric current passes through the true contact area between contacting solids subjected to a difference of the electric potential. However, due to the presence of non-conducting oxides on the surface the effective conducting area is, in general, smaller than the true contact area. (6) Adhesion between solids is also considerably controlled by roughness and the effective contact area. On the other hand, the adhesion exploited by insects and small animals such as gecko [@scherge2001biological; @huber2007influence] is more related to the hierarchical structure of their limbs. (7) Fluid flow between contacting (or slightly separated) surfaces in sealing applications is controlled by their roughnesses, the out-of-contact areas, forming an opening field, and enabling the fluid to pass through the interface. Under certain critical load no possible paths exist, which corresponds to the percolation point. (8) Frictional sliding on contacting asperities is a dissipative process, the produced heat partly goes into the rubbing solids and may results in huge thermal gradients and very high localized temperatures at sliding asperities [@bowden_tabor; @rice2006heating]. These extreme thermo-mechanical loads can induce phase transformations and changes in material microstructures. (9) Local wear and initiation of fatigue cracks in fretting conditions [@proudhon2005ijf; @kubiak2011ti] strongly depend on the local stress state, which is intimately linked with the surface roughness. (10) Produced by adhesive wear [@sasada1984w], debris particles are necessarily related to the true contact area as discussed in (3); geometrical characteristics of surfaces are less important for the debris formed in abrasive wear [@hokkirigawa1988ti]. The presence of debris (or mobile third body) can significantly affect the true contact area and thus modify the kinetics of debris production. (11) Effect of the true contact area on the boundary lubrication is straightforward, concerning the hydrodynamic lubrication regimes and mixed lubrication the opening (gap) field plays an essential role and determines the film thickness, its possible breakage and thus determines the resulting friction and wear. In all the describing processes and phenomena the mechanical interaction between rough surfaces plays the key role. Accurate prediction of this mechanical interaction would allow constructing reliable and physically sound models for all the aforementioned phenomena. Mechanical models of rough contact can be divided into two main classes: deterministic and statistical models. The former models rely on direct numerical simulations of contact [@pei2005jmps; @campana2006practical; @anciaux2009ijnme; @yastrebov2011rough; @putignano2012ijss; @afferrante2012w; @bemporad2015optimization; @muser2017meeting] and take as an input the topographical data of contacting surfaces and their material properties in order to predict integral quantities such as interface stiffness, evolution and growth of true contact area and of the mean gap under increasing pressure, etc. In addition, such models allow to study the shape and spatial distribution of contacting clusters, local stress states, critical zones, etc. These models are especially important for situations when the contact area is small or when the roughness of contacting solids cannot be seen as homogeneous at the scale of contact area, for example, the rough surface can have critical defects, rare high asperities or rare third-body particles. On the contrary, the statistical models [@greenwood1966prcl; @bush1975w; @persson2001jcp; @carbone2008jmps] require only a few parameters of surface roughness (for example, root-mean squared height or surface gradient, probability density, fractal dimension) and are appropriate for situations of homogeneous roughness. These models can yield results in terms of statistical distributions (for example, pressure and tangential tractions) and mean integral values (for example, mean spacing between contacting asperities, integral contact area). Such models are, in general, easier to use than deterministic but on the other hand they are inherently less precise since cannot take into account finite size of the contact area and cannot accurately resolve complex mechanical interaction between contacting spots.

In this paper we consider a contact problem between two non-conformal solids with superposed roughness. This generalization of Hertzian contact is relevant to almost all aforementioned applications since in real engineering and natural systems, the contact can be usually reduced to a contact between two smooth surfaces, whose shape is described by a two-form. Here however, we limit our attention to the contact of solids of revolution. To address this problem which was already studied using different models in different regimes [@greenwood1967elastic; @goryacheva2006mechanics; @goryacheva2013contact; @pastewka2016contact], we choose to use a deterministic multi-asperity approach, rather similar to ones elaborated in [@afferrante2012interacting; @afferrante2018elastic], which includes long-range elastic interaction between asperities. The main objective is to demonstrate capacities of the model in predicting particular features of this contact interaction for a given (deterministic) roughness. On the other hand we will explore its ability to derive statistical results, which, however, can be also done using Greenwood and Tripp’s model [@greenwood1967elastic]. The deterministic multi-asperity model, compared to statistical multi-asperity models [@greenwood1966prcl; @bush1975w; @carbone2008jmps], takes into account long-range interactions; and it enables to handle accurately very light loads, which are hardly reachable for conventional direct numerical simulations using finite element [@pei2005finite; @yastrebov2011cras] or modifications of boundary element methods [@stanley1997; @polonsky1999numerical; @campana2006practical; @bemporad2015optimization].

The paper is organized as follows. In Section \[sec:1\] we define the problem, recall the roughness generation method and explain the method to extract asperities’ data. In Section \[sec:2\], the numerical model to solve contact problems is presented. The results are presented in Section \[sec:3\], namely (1) the analysis of the true contact area evolution, the two relevant regimes and transition between them; (2) the study of contact pressure and contact area as well as their spatial and statistical distributions. We conclude by a short discussion of main results in Section \[sec:4\].

\[sec:1\]Problem set-up and roughness
=====================================

Set-up
------

Using the deterministic multi-asperity model with long-range interaction we will consider a problem of contact between locally axisymmetric parabolic solids \[see Fig. \[fig:setup\](a)\] of curvature radii $R_1$ and $R_2$, the solids are assumed to be elastic, homogeneous and isotropic with elastic properties $E_1,\nu_1$ and $E_2,\nu_2$, respectively, where $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio, respectively. Schematically this problem can be presented as a contact between two spheres[^1], see Fig. \[fig:setup\]. We assume that the contacting solids possess rough surfaces, which, in the close vicinity of contact zone, can be described by analytical functions $z_1(x,y)$ and $z_2(x,y)$. Note that this roughness in case of curved surfaces shall be superimposed on the macroscopic shapes, which can be assumed to be $z_{si}=(x^2+y^2)/(2R_i)$. We assume that we can neglect the tangential (in-contact-plane) displacement of surface points, i.e. the effective roughness is given by $z^*(x,y) = z_1(x,y)-z_2(x,y)+z_0$, where $z_0$ is a constant. Under frictionless elastic contact between rough surfaces with relatively small slopes, this problem can be mapped on a problem between a smooth rigid sphere with effective radius $R^* = R_1R_2/(R_1+R_2)$ and a rough elastic half-space with effective elastic modulus $E^*$ [@johnson1985b; @barber2003bounds]: $$E^* = \frac{E_1E_2}{(1-\nu_2^2)E_1 + (1-\nu_1^2)E_2}.$$ In case of smooth surfaces, we recall [@johnson1985b] below the analytical solutions for this problem under the action of a squeezing force $N$. The contact radius $a$, the approach between distance points $\delta$, and the pressure distribution $p(x,y)$ are given by the following expressions: $$a = \left(\frac{3NR^*}{4E^*}\right)^{1/3},\quad \delta = \frac{a^2}{R^*} = \left(\frac{9N^2}{16R^*E^{*2}}\right)^{1/3},\quad p(x,y) = \frac{3P}{2\pi a^2} \sqrt{1-\frac{x^2+y^2}{a^2}}.$$

![\[fig:setup\]Contact between two elastic spheres with rough surfaces and the equivalent problem of indentation of an effective elastic and rough half-space by a rigid smooth sphere of effective radius.](setup){width="100.00000%"}

Roughness
---------

To model the surface roughness we use a white-noise filtering technique [@hu1992ijmtm; @yastrebov2017jmps]. This technique enables us to control accurately the spectrum of the surface and also to preserve the stochastic aspect of real roughness. Namely, we control (1) the Hurst exponent $H$, which determines the power-law decay (at least for short wavelengths $\lambda$) of spectral energy $\Phi \sim \lambda^{2(1+H)}$, and also (2) the lower and (3) the upper cut-off wavelengths $\lambda_s$ and $\lambda_l$, respectively, where lower indices $s$ and $l$ stand for “short” and “long”. At the same time the spectrum can be made bi-fractal: a power-law increase with the wavelength in the interval $\lambda\in[\lambda_s,\lambda_l]$ can be followed by a plateau for longer wavelengths $\lambda\in(\lambda_l,L]$, where $L$ is the sample length.

Based on this definition, we can readily conclude that for loads implying the value of the macroscopic contact radius ${a_{\mbox{\tiny macro}}}$ which is comparable with the shortest roughness wavelength ${a_{\mbox{\tiny macro}}}\sim\lambda_s$, the contact will be formed by a few contact spots, whereas for ${a_{\mbox{\tiny macro}}}\gg \lambda_s$ it is expected that the contact contains a statistically meaningful number of spots. Therefore, statistical models could be accurately applied only for the later case. At the same time, it is known that the longest wavelength (of the power-law spectral decay) $\lambda_l$ determines the dispersion of results [@yastrebov2012pre; @pastewka2013pre; @yastrebov2015ijss; @rey2017hal]: therefore, to have a spatially homogeneous and statistically reproducible contact-spot distribution for all possible roughness realizations, the contact radius has to fulfil the following condition ${a_{\mbox{\tiny macro}}}\gg \lambda_l$ as well. However, this last condition should be fulfilled only for generated rough surfaces possessing discrete spectra.

Another important consideration is that in multi-asperity models it is assumed that local contact spots are located far from each other. A softer restriction would be the assumption that the contact zones associated with different asperities do not intersect. This strong limitation was partly overcame in [@afferrante2012w] by assuming that coalescing contact spots (non-convex and possibly non-simply connected) form a new effective one of a circular shape. However, this improvement requires to switch from statistical to a deterministic model. Moreover, separate asperities located on a smooth surface are inevitably separated by valleys and their merge is associated with a formation of junctions over saddle points [@yastrebov2014tl; @dapp2015epl], which cannot be described in the framework of Hertzian theory which assumes the parabolic shape of contacting solids. Rigorously in the framework of statistical model we shall assume that at least the contact radius associated with asperities does not overpass the quarter[^2] of the average distance between asperities $\langle d\rangle$ (see \[app\] for details), i.e. $${a_{\mbox{\tiny micro}}}< \frac{\langle d\rangle}{4} = \frac{1}{4\sqrt D} \approx \lambda_l \sqrt{\frac{\sqrt3 }{8  \pi } \frac{(2-H)}{(1-H)}\frac{(\zeta^{2-2H} - 1)}{\zeta^{4-2H}}}   \label{amicro_limit}$$ where $D$ is the density of asperities, and the last approximate equality is valid for surfaces with a rich spectral content, i.e. for $\zeta = \lambda_l/\lambda_s \gg 1$.

In most statistical models, both multi-asperity and Persson’s model, the rough surface is assumed to be Gaussian, i.e. the heights have an infinite support. A notable exception is the study of asperity-height distribution on a semi-infinite support [@ciavarella2017asme]. Therefore, for any reasonable approaching (here, the distance between mean planes of rough surfaces expressed in roughness-dispersion units), a small but non-zero portion of the surface will experience an infinite displacement, which would readily result in an apparent violation of criterion . However, since in statistical models the notion of asperities is hidden behind the height probability density function, it seems reasonable to find the mean local asperity radius by assuming that the contact area ${A_{\mbox{\tiny asp}}}$ predicted by a multi-asperity model over nominal area $A_0$ is formed by $N$ identical contact spots $\pi{a_{\mbox{\tiny micro}}}^2$, where the number of spots $N$ is assumed to be known through the asperity density $D$ as $N = DA_0$, therefore the individual contact-spot radius can be estimated as $${a_{\mbox{\tiny micro}}}\approx \sqrt{\frac{{A_{\mbox{\tiny asp}}}}{\pi D A_0}},$$ where from using , a simple limit for the validity of multi-asperity models could be established: $$\frac{{A_{\mbox{\tiny asp}}}}{A_0} < \frac{\pi}{16} \approx 19.6\%.$$ This limit can be used locally (point-wise) as the upper bound for multi-asperity models, i.e. the local area fraction predicted by multi-asperity model cannot overpass $19.6$ %. Of course, because of non-uniform distribution of asperities and contact spots, the real limit can be considerably smaller (see also comparison between deterministic and multi-asperity models [@yastrebov2017jmps]).

Asperities
----------

In vicinity of every summit, a smooth surface can be approximated by a quadratic form of local coordinates: $$z(x,y) = z_0 + a x^2 + b y^2 + 2 c xy = z_0 + {\left[\begin{array}{c}}x\\y{\end{array}\right]}^\intercal {\left[\begin{array}{cc}}a \;\;& c \\ c \;\; & b{\end{array}\right]}{\left[\begin{array}{c}}x\\y{\end{array}\right]}= z_0 + {\boldsymbol{X}}^\intercal {\boldsymbol{A}} {\boldsymbol{X}}
$$ with $$a = \frac12{\frac{\partial^2{z}}{\partial{x}^2}}, \quad b = \frac12{\frac{\partial^2{z}}{\partial{y}^2}}, \quad c = \frac12\frac{\partial^2 z}{\partial x\partial y}.$$ To determine the surface principal curvatures and their orientation, we seek for such rotation $\varphi$ of $\{x,y\}$ coordinates $${\boldsymbol{X}} = {\boldsymbol{Q}} {\boldsymbol{X}}',\quad {\boldsymbol{X}} =  {\left[\begin{array}{c}}x\\y{\end{array}\right]},\quad {\boldsymbol{X}}' =  {\left[\begin{array}{c}}x'\\y'{\end{array}\right]},\quad {\boldsymbol{Q}} =  {\left[\begin{array}{cc}}\cos(\varphi) & \sin(\varphi)\\ -\sin(\varphi) & \cos(\varphi) {\end{array}\right]}$$ that the resulting form is diagonal: $$z(x,y) = z_0 + {\boldsymbol{X}}'^\intercal {\boldsymbol{Q}}^\intercal {\boldsymbol{A}} {\boldsymbol{Q}} {\boldsymbol{X}}', \quad {\boldsymbol{Q}}^\intercal {\boldsymbol{A}} {\boldsymbol{Q}} = {\left[\begin{array}{cc}}\lambda_1 \;\; & 0 \\ 0\;\; & \lambda_2 {\end{array}\right]}.$$ The angle $\varphi$ that ensures the diagonality should verify the following equation: $$(a - b) \cos(\varphi)\sin(\varphi) + c (\cos^2(\varphi) - \sin^2(\varphi))  = 0,$$ which after dividing by $(c\cos^2(\varphi))$ gives a quadratic equation for $\tan(\varphi)$: $$\tan^2(\varphi) - \frac{a-b}{c}\tan(\varphi)  - 1 = 0,$$ therefore the two solutions are given by : $$\varphi_{i} = \mathrm{arctan}\left(\frac{a-b}{2c} \pm \sqrt{\left(\frac{a-b}{2c}\right)^2+1}\,\right),
 \label{eq:angle}$$ which correspond to eigenvalues $$\lambda_{i} = \frac12\left(a+b \pm \sqrt{(a-b)^2 + 4c^2}\right),$$ therefore the principal curvatures and principal radii are given by: $$\kappa_{i} = \left(a+b \pm \sqrt{(a-b)^2 + 4c^2}\right),\quad R_{i} = 1\left/\left(a+b \pm \sqrt{(a-b)^2 + 4c^2}\right)\right.
  \label{eq:curvatures}$$

Rough surface measurements are inevitably of discrete nature, therefore the question of asperity identification is relevant [@paggi2010w]. Numerically generated rough surfaces [@hu1992ijmtm; @meakin1998fractals] (as soon as they are not given by analytical functions) are also discrete, i.e. the surface roughness is determined by a matrix of heights $z_{ij}$, where lower indices $i,j$ determine the in-plane position as $z_{ij} = z(i\Delta x, j\Delta y)$, where $\Delta x,\Delta y$ are the sampling intervals in $x,y$ directions, respectively. For such a representation, every potential summit $z'_{ij}$ is identified if the nearest and next-nearest (diagonal) neighbour points are lower: $z'_{ij} > z\left[(k+i)\Delta x,(l+j)\Delta y\right]$, where $k,l$ may take values $\{-1,0,1\}$ except $\{k,l\}=\{0,0\}$. Using the finite differences, for all generalized summits the following values of the components of the fundamental matrix ${\boldsymbol{A}}$ can be estimated as: $$a = \frac12{\frac{\partial ^2\!z}{\partial x^2}} \approx \frac{z_{i+1j}+z_{i-1j}-2z_{ij}}{2\Delta x^2};\quad b=\frac12{\frac{\partial ^2\!z}{\partial y^2}} \approx\frac{z_{ij+1}+z_{ij-1}-2z_{ij}}{2\Delta y^2}$$ $$c = \frac12{\frac{\partial ^2\!z}{\partial x\partial y}} \approx \frac{z_{i+1j+1}+z_{i-1j-1}-z_{i+1j-1}-z_{i-1j+1}}{8\Delta x \Delta y}.$$ The principal curvatures and radii then can be found from . If both curvatures are negative the generalized summit is the real summit. If the principal curvatures have a different sign $\kappa_1\kappa_2 < 0$ this “summit” is a saddle point masked by discretization. The algebraic $\kappa_a$ and the geometric $\kappa_g$ mean curvatures are defined as follows $$\kappa_a = \frac12(\kappa_1+\kappa_2) = \frac12(a+b)$$ $$\kappa_g = \sqrt{\kappa_1\kappa_2} = 2\sqrt{ab - c^2}$$

![\[fig:asperities\]Generated rough surfaces (upper panel), detected elliptic asperities (middle panel) and approximate circular asperities with geometric mean curvature. For all surfaces, the Hurst exponent $H=0.8$, $L/\Delta x=2048$, lower cut-off $L/\lambda_l=4$, and root mean squared height $\sigma/L=0.025$: (a) $L/\lambda_s=16$, (b) $L/\lambda_s=32$, (c) $L/\lambda_s=64$. The surfaces are generated from the same white noise, therefore they are visually similar. The colour of asperities is selected according to the height of the summit.](asperities){width="100.00000%"}

In Fig. \[fig:asperities\] we show generated rough surfaces with different spectral content and the detected asperities approximated by ellipses with relevant orientation and circles with the geometric mean curvature. Of course, the ellipses and circles give only an approximate representation of the roughness: anyway, the Hertzian theory is valid for arbitrary elliptic paraboloids, which can be approximated by ellipsoids near their summits. The size of ellipses and circles corresponds to $L=1$ (l.u.[^3]) and root mean squared height (RMS height) $\sqrt{\langle (z-\bar z)^2\rangle}/L = \sigma/L =  0.025$. Obviously, a decrease of the RMS height would result in an increase of asperities’ radii. It was shown by Greenwood [@greenwood2006w], that according to Nayak’s theory [@nayak1971tasme], the asperities are only “mildly” elliptic, and thus, the geometric mean curvature can be successfully used in multi-asperity contact models [@greenwood1985elliptic] instead of treatment elliptic contacts [@bush1975w; @johnson1985b], without any significant loss in accuracy. However, as seen by “naked eye” from Fig. \[fig:asperities\], a considerable population of asperities has a pronounced ellipticity. Nevertheless, the mechanical analysis can be successfully carried out using simple geometric mean curvature and associated correction factors even for strongly elliptic asperities, see Table 2 in [@greenwood2006w].

It seems that a good test for the surface generator would be a verification of agreement between asperities’ statistics detected on the surface and Greenwood’s analysis [@greenwood2006w], who obtained the following equation for the joint probability density of the ensemble of asperities’ principal curvatures[^4]: $$P(\kappa'_1,\kappa'_2) = \frac{27}{16\sqrt\pi}\kappa_1'\kappa_2'|\kappa_1'-\kappa_2'|\exp\left[-\frac{3}{16}\left(3\kappa_1^{\prime 2} + 3\kappa_2^{\prime 2} - 2\kappa'_1\kappa'_2\right)\right],
  \label{eq:Greenwood}$$ where $\kappa_i' = \kappa_i/\sqrt{m_4}$ represent normalized curvatures. This comparison is shown in Fig. \[fig:asp\_stat\]. Numerically, the statistics is performed over an ensemble of asperities detected on a single realization of a surface with a rich spectrum $L/\lambda_l = 4$, $L/\lambda_s=2048$, $L/\Delta x = 8192$ and $H=0.8$; in total this surface contains $\approx 580\,000$ asperities. In accordance with Greenwood’s findings, the most probable asperities possess the following ratio of principal curvatures $\kappa_1/\kappa_2 \approx 2.214$; the probability to find an asperity close-to-spherical is marginal. The overall mean geometrical curvature for the ensemble of asperities can be obtained from  by integration $$\bar\kappa'_g = \int\limits_0^\infty\int\limits_0^\infty \sqrt{\kappa_1'\kappa_2'} P(\kappa'_1,\kappa'_2) \,d\kappa_1'd\kappa_2' \approx 1.356
  \label{eq:mean_kg}$$

![\[fig:asp\_stat\]Comparison of the joint probability density computed numerically (on the left) and analytically, following Greenwood’s equation ; principal curvatures $\kappa_1,\kappa_2$ are normalized by the average algebraic curvature $\bar\kappa_a$.](Greenwood_statistics){width="100.00000%"}

\[sec:2\]Numerical model
========================

The numerical model is based on a deterministic multi-asperity Hertzian contact with included long range elastic interaction between asperities. A set of asperities $\mathcal I$ with in-plane coordinates $x_i,y_i$, summit height $z_i$ and curvature radius $r_i$ is distributed on the top of an elastic half-space. The asperity parameters come directly from the generated rough surfaces. This set of elastic asperities comes in contact with a rigid sphere of radius $R$. The sphere’s tip is located at $\{X_0,Y_0,Z_0\}$. Then the penetration between the sphere and $i$-th asperity is determined by: $$g_i = z_i  - \Delta z_i +\sqrt{R^2-(X_0-x_i)^2-(Y_0 -y_i)^2} - R - Z_0 = g_{i0} - \Delta z_i ,
  \label{eq:pen}$$ where $\Delta z_i$ is the vertical displacement of the asperity induced by contact forces acting on other asperities of the surface, and $g_{i0} = z_i  +\sqrt{R^2-(X_0-x_i)^2-(Y_0 -y_i)^2} - R - Z_0$. If the penetration is positive, Hertzian contact theory can be used to determine the resulting force $f_i$ and contact radius $a_i$ as: $$f_i = c_i\langle g_i\rangle^{3/2}, \quad a_i = \sqrt{r_i  \langle g_i\rangle}, \quad\mbox{ with } \quad c_i = \frac43 E^* \sqrt{r_i}   \label{eq:force}$$ where $\langle g_i\rangle = \max\{0,g_i\}$ denotes Macaulay brackets, therefore the forces are induced only if the penetration is positive. We assume that the half-space’s surface deforms under action of contact forces following Boussinesq solution, thus the vertical displacement due to a single concentrated force $f$ located at $x_f,y_f$ is given by: $$\Delta z_f = \frac{f}{\pi E^*\sqrt{(x-x_f)^2+(y-y_f)^2}}.
  \label{eq:dz}$$ This coupling between additional displacement and forces makes the system of equations - strongly coupled. Then, following the superposition principle, at a given asperity $i$, its total vertical shift is given by: $$\Delta z_i = \frac{1}{\pi E^*}\sum\limits_{j \in \mathcal I\setminus i} \frac{f_j}{d_{ij}} =  {\boldsymbol{B}} {\boldsymbol{F}},$$ where $d_{ij} = \sqrt{(x_i-x_j)^2+(y_i-y_j)^2}$ and the symmetric interaction matrix ${\boldsymbol{B}}$ (containing the inverse distances) and the vector of forces ${\boldsymbol{F}}$ are given by: $${\boldsymbol{B}} = 
  \frac{1}{\pi E^*}\left[\begin{array}{ccccc}\displaystyle        0 & d^{-1}_{12} & d_{13}^{-1} & \quad\dots & d_{1N}^{-1}\\[3pt]
        d^{-1}_{12} & 0 & d_{23}^{-1} & \quad\dots & d_{2N}^{-1}\\[3pt]
        d^{-1}_{13} &  d_{23}^{-1} & 0 & \quad\dots & d_{2N}^{-1}\\[3pt]
        \vdots & \vdots  & \vdots  & \ddots & \vdots \\[3pt]
        d^{-1}_{1N} \quad &  d_{N2}^{-1}\quad & d_{N3}^{-1} \quad& \quad\dots & 0\\[3pt]
           \end{array}
           \right], \quad {\boldsymbol{F}} = \left[\begin{array}{c}
                                          f_1\\[3pt]f_2\\[3pt]\\[3pt]\vdots\\[3pt] f_N
                                         \end{array}\right]$$ Therefore from  we obtain a non-linear system of equations: $${\boldsymbol{F}} - {\boldsymbol{C}} \langle {\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}}\rangle^{3/2} = 0,
 \label{eq:sys}$$ where ${\boldsymbol{C}} = \mathrm{diag}[c_1 \; c_2 \; \dots \; c_N]$ is the diagonal matrix, and the vector ${\boldsymbol{G}}_0 = [g_{10} \; g_{20} \; \dots \; g_{N0}]^\intercal$ contains the local penetrations in absence of elastic interactions. In used notations we also assume that Macaulay brackets applied to a vector and raising to $n$-th power imply the following operation with vector components: $$\langle {\boldsymbol{V}} \rangle = [\langle {\boldsymbol{v}}_1 \rangle\; \langle {\boldsymbol{v}}_2 \rangle\; \dots \; \langle {\boldsymbol{v}}_N \rangle\;]^\intercal,\quad {\boldsymbol{V}}^n = [v_1^n \; v_2^n \; \dots \; v_N^n]^\intercal.$$ Iterative explicit resolution of this problem (i.e. when to compute the configuration on the current load step, forces from the previous load steps are used to estimate the long-range interaction) in most cases leads to big errors. On the other hand, in implicit integration, the system converges badly if Newton-Raphson method is used, because of the discontinuous nature of the function $\langle x\rangle$ and non-coercive properties. Therefore we regularize it as follows $$m(x,\delta) =
  \left\{ \begin{array}{ll}
          \displaystyle\frac{\delta^2}{2\delta-x} ,&\mbox{ if } x < \delta\\
          x,&\mbox{ otherwise},
            \end{array}\right.$$ If $\delta \to 0$ we recover the classical Macaulay brackets: $m(x,\delta) \to \langle x\rangle$. This approximation creates small repulsive forces between the sphere and asperities, but they decay quite rapidly (as $z^{-3/2}$), the parameter $\delta$ reduces to zero during the convergence loop enabling to ensure a better convergence for Newton-Raphson method, initial value of $\delta$ is chosen to be $\sigma/100$. With this smoothing the system of equations  reduces to $${\boldsymbol{F}} - {\boldsymbol{C}} m( {\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}},\delta)^{3/2} = 0,
 \label{eq:sys_approx}$$ whose solution tends to the solution of  when $\delta\to0$. The initial guess for the force vector ${\boldsymbol{F}}_0$ is computed from a model without interaction (${\boldsymbol{B}} = 0$), then from Eq. : $${\boldsymbol{F}}_0 = {\boldsymbol{C}} \langle {\boldsymbol{G}}_0 \rangle^{3/2}.$$ A classical choice to take the initial guess from previous converged load step is not used since the possible set of active asperities changes from step to step (see next Section) and even if not, such initial guess often leads to divergence of the computational step. From the initial guess, we obtain the increment of $\Delta {\boldsymbol{F}}$ as: $$\Delta{\boldsymbol{F}} = -{\boldsymbol{K}}^{-1}{\boldsymbol{R}},$$ where the free term ${\boldsymbol{R}}$ and the tangent matrix ${\boldsymbol{K}}$ are given by: $${\boldsymbol{R}} = {\boldsymbol{F}}_0 - {\boldsymbol{C}} m( {\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}}_0,\delta)^{3/2},$$ $${\boldsymbol{K}} = {\boldsymbol{I}} - \frac32 {\boldsymbol{C}} m({\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}}_0,\delta)^{\frac12}\left.{\frac{\partial m({\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}},\delta)}{\partial {\boldsymbol{F}}}}\right|_{{\boldsymbol{F}} = {\boldsymbol{F}}_0},$$ where ${\boldsymbol{I}}$ is identity matrix. After every iteration, according to the classical Newton’s method, the vector of forces is updated as ${\boldsymbol{F}}^i = {\boldsymbol{F}}^{i-1}+\Delta{\boldsymbol{F}}$ and all negative forces are removed as ${\boldsymbol{F}}^i = \langle {\boldsymbol{F}}^i\rangle$. Regardless the fact that function $m(x,\delta)$ is smooth and well behaved, the convergence of the Newton method is still slow and non-monotonous, but it is rather robust. The convergence is assumed to be reached when $$\|{\boldsymbol{\Delta}} {\boldsymbol{F}}^i \| / \|{\boldsymbol{F}}^i \| \le \varepsilon,$$ the tolerance in all simulations is chosen to be $\varepsilon = 0.1$%.

Computational set-up
--------------------

We assume that the spherical indenter of radius $R$ is rigid and that its tip is centred at $x=X_0$, $y=Y_0$ of in-plane coordinates. The substrate and asperities are made of gold with the following average elastic properties: the Poisson’s ratio $\nu=0.42$ and the Young’s modulus is $E=78$ GPa, giving the effective modulus $E^* \approx 94.71$ GPa. The size of the square region on elastic half-space with identified asperities is set to be $1$ mm $\times$ $1$ mm, i.e. the maximal contact radius is limited to $0.5$ mm. The indentation simulations are carried out over a set of deformable asperities located on a deformable half-space; very small indentation depths are used. Naturally, no periodic boundary conditions are used. The sphere is approached in 100 load steps from the separation[^5] $u_z = 4\sigma$ down to $u_z = \sigma$ in 100 load steps. In order to accelerate the simulations, only those asperities are considered which (1) “penetrate” into the sphere in the undeformed state, i.e. $z_i > u_z + R-\sqrt{R^2 - (X_0-x_i)^2 - (Y_0-y_i)^2}$, and (2) which are located in the circle $\sqrt{(x-X_0)^2 + (y-Y_0)^2} < r_0$, where $r_0 = 1.5\sqrt{2(4\sigma - u_z)R}$, which is 1.5 times bigger than the radius of geometrical overlap of a sphere, whose tip is penetrated under a plane by distance $(4\sigma-u_z)$.

After convergence, at every step we record the total contact force $N$ and total contact area $A$. Moreover, at every load step, a coarse grained map of the contact pressure and the contact area is saved. The map represents a regular grid $N\times N$ created over the indentation area $L\times L$ with cells of size $\Delta x \times \Delta x$. In each cell $\{i,j\}$, $i,j \in [0,N-1]$ all contacting asperities $\mathcal I_{ij} \subset \mathcal I$ contribute to the coarse-grained contact area fraction $A_{ij}$, and to the coarse-grained contact pressure $p_{ij}$, where $$A'_{ij} = \frac{1}{\Delta x^2}\sum\limits_{k\in\mathcal I_{ij}} \pi a_k^2,\quad p_{ij} = \frac{1}{\Delta x^2} \sum \limits_{k\in\mathcal I_{ij}} F_{k},
  \label{eq:coarse}$$ where $a_k$ is computed using  and $\mathcal I_{ij}$ is a set of asperities with in-plane coordinates located in the given cell.

\[sec:3\]Results
================

Two regimes of the contact-area growth
--------------------------------------

From the physical evidence, under most conditions and for most materials, the frictional force (static friction) increases practically linearly with the normal load independently on the nominal contact area, which is reflected by the Amontons-Coulomb’s friction law. Accordingly to the adhesive theory of contact [@bowden_tabor], the frictional force is proportional to the true contact area. Therefore, the linear growth of the contact area could be considered as a necessary condition of the model validity. The model of Greenwood and Tripp [@greenwood1967elastic], which represents the fully analytical treatment of the elastic indentation problem of a rough surface including elastic interaction between contacting asperities in the statistical sense , predicts the linearity between the force and the contact area, as well as the recent full scale numerical treatment of this problem, which was carried out by Pastewka & Robbins [@pastewka2016contact]. However, as remarked in [@greenwood1970contact; @pastewka2016contact], at very small loads/large separations a single-asperity contact regime exists, which results in Hertzian scaling for elastic asperities[^6].

Expectedly, our numerical study with the deterministic multi-scale model identifies these two deformation regimes (see Fig. \[fig:area\]). The data were averaged over 120 simulations per each combination of root mean squared roughness $\sigma=\{0.001,0.01,0.1\}$ and different upper cut-offs $L/\lambda_s=\{256,512\}$, the lower cut-off was kept constant $L/\lambda_l = 24$. The vertical displacement of the indenter’s tip was changed from $Z_0 = 4\sigma$ to $\sigma$ in 100 load steps, the indenter’s radius is fixed to be $R=2.5$ mm. The error bars correspond to a single standard deviation; the fact that error bars are missing in very low pressure-regime is explained by the fact that only in a single realization the contact occurs at such small loads.

The first deformation regime corresponds to the deformation of a single Hertzian asperity, it occurs at very light loads and the contact-area grows as $A_I \sim N^{2/3}$. The second regime is characterized by an approximately linear evolution of the contact area $A_{II}\sim N$ and it corresponds to multi-asperity statistical contact. The first regime then fully depends on the curvature of the first contacting asperity, which can be assumed to be the average geometrical mean curvature over all asperities[^7] $r = \bar r_g \approx 0.737/\sqrt{ m_4}$ (the last approximate equality follows from ): $$A_I = \pi \left(\frac{3\bar r_g}{4E^*}\right)^{\!\!\frac 2 3} N^{^{\frac 2 3}} \approx \pi \left(\frac{0.55275}{\sqrt{ m_4}E^*}\right)^{\!\!\frac 2 3} N^{^{\frac 2 3}}.
\label{eq:AI}$$ In the second regime, at moderate loads, the contact area should grow accordingly to classical multi-asperity theories [@greenwood2006w; @bush1975w; @carbone2008jmps], Persson’s theory [@persson2001jcp] and more recent full-scale numerical simulations [@hyun2004pre; @pei2005jmps; @putignano2012jmps; @afferrante2012w] with the scaling depending on the root mean squared roughness gradient $\langle |\nabla z|^2\rangle$, which for Gaussian isotropic surface can be estimated through the second spectral moments $\langle |\nabla z|^2\rangle = \sqrt{2m_2}$ (see \[app\]). Therefore, in the second regime the contact area growth is determined by the following equation: $$A_{II} = \frac{\beta\kappa}{\sqrt{2 m_2}E^*} N,
\label{eq:AII}$$ where the constant $\beta \approx 0.2$ was identified from our simulations to fit approximately the data. The factor $\kappa=\sqrt{2\pi}$ appears from statistical multi-asperity models [@bush1975w; @greenwood2006w; @carbone2008jmps] and was shown to be very accurate for predicting the contact area evolution at small loads for nominally flat surfaces [@yastrebov2017jmps]. For higher loads, both multi-asperity models [@carbone2008jmps; @greenwood2006w] and numerical simulations [@paggi2010w; @yastrebov2017jmps] demonstrate a slightly non-linear area growth which depends on Nayak’s parameter $\alpha = m_0m_4/m_2^2$: greater $\alpha$ results in smaller contact area. Note that in our simulations we remain far from the third deformation regime, purely Hertzian one, which was identified in [@pastewka2016contact].

![\[fig:area\]Evolution of the contact area for different root mean squared roughness $\sigma=\{0.001,0.01,0.1\}$ and different upper cut-offs $L/\lambda_s=\{256,512\}$; to guide the eye, the $2/3$-power law and a linear evolution with relevant factors (Eqs. ,) are also plotted.](Force_area_average_non_norm){width="100.00000%"}

Transition between the regimes
------------------------------

The transition force $N^*$ between the two regimes can be found by assuming equality of areas $A_I(N^*) = A_{II}(N^*)$, which results in the following expression: $$N^* = 0.55275^2 \frac{\pi^3 (2  m_2)^{3/2} E^*}{\beta^3\kappa^3  m_4} \approx 212.66  \frac{  m_2^{3/2} E^*}{ m_4},
  \label{eq:trans}$$ where the numerical factor $212.66 \approx 0.55275^2 (\pi\sqrt2/\kappa \beta)^3$ with $\kappa = \sqrt{2\pi}$ and $\beta = 0.2$. Note that in [@pastewka2016contact], a different expression was obtained[^8]: $N^*_{\mathrm{PR}} = (9\pi^3/4m_4) (\sqrt{2 m_2}/\kappa)^3E^*$, which after evaluating all numerical quantities and assuming $\kappa = \sqrt{2\pi}$ gives $N^*_{\mathrm{PR}} \approx 12.5287 m_2^{3/2}  E^*/m_4$ which is considerably smaller than our value . This is because of slight overestimation made by the authors [@pastewka2016contact] of the asperity radius $\bar r_g = 2/\sqrt{m_4}$ and of non-taking into account the scaling factor[^9] $\beta$ which was required in our case to predict the contact area growth in the second regime. On the other hand, Pastewka & Robbins’ expression in its final form lacks the fourth spectral moment $m_4$, which in their case was expressed through the short wavelength cut-off, Hurst exponent and the second spectral moment. Returning to our results, we show that the lower cut-off $\lambda_l$ affects slightly the numerical values, but keeps the transition interval unchanged. The normalized data plotted in Fig. \[fig:area\_norm\] include the results found for different indenter’s radius $R=\{2.5;\;5\}$ mm, and different lower cut-off wavelength $\lambda_l/L = \{24;\;48\}$, which seemingly do not change the transition point when normalized.

Scaling
-------

Since all spectral moments scale proportionally to the roughness variance $m_4 \sim m_2 \sim m_0 = \sigma^2$ then the contact areas in the first  and the second  regimes scale as follows: $$A_I(N,\sigma) = A_I(N/\sigma) \sim \left(\frac{N}{\sigma}\right)^{2/3},\quad A_{II}(N,\sigma) =  A_{II}(N/\sigma) \sim \frac N \sigma.$$ The scaled results and theoretical predictions are plotted in Fig. \[fig:area\_norm\]. Note that the scaling by the RMS gradient is used which is equivalent. In the geometrical model, scaling the RMS roughness scales proportionally the RMS gradient and curvature, however, it should be bared in mind that in real life, changing RMS roughness (removing higher asperities), in general, does not result in changing the RMS gradient which depends on small scale roughness, therefore the latter was used to scale the results. According to this scaling, the transition load between the two regimes should be proportional to the RMS height $N^* \sim \sigma$.

![\[fig:area\_norm\]Evolution of the contact area as a function of the force normalized by the RMS gradient, the data are obtained for different root mean squared roughness $\sigma=\{0.001,0.01,0.1\}$ and different upper cut-offs $L/\lambda_s=\{256,512\}$; if it is not stated differently $L/\lambda_l = 24$ and $R=2.5$ mm; to guide the eye, the $2/3$-power law and a linear evolution with relevant factors (Eqs. ,) are also plotted.](Average_force_area_normalized){width="100.00000%"}

Maps of contact area/pressure
-----------------------------

The linear regime  seemingly follows directly from classical theories. If one assumes that a smooth Hertzian solution can be used at macro scale, and that at the micro scale the contact area is proportional to the local contact pressure, then the total contact area is proportional to the normal force. Let the local contact area $dA$ be given by $dA(x,y) = \kappa p(x,y)dA_0/(\sqrt{2m_2}E^*)$ and the pressure is distributed according to Hertz theory $p(x,y) = p_0\sqrt{1-(x^2+y^2)/a^2}$, with $p_0 = 3N/(2\pi a_0^2)$, where $N$ is the squeezing force. Since the integral of pressure over the contact area is simply the resulting force $N$, then the contact area is given by: $$A_{\mathrm{th}} = \frac{\kappa}{\sqrt{2m_2}E^*}N .
  \label{eq:Ath}$$ However, the contact area for indentation of the rough contact [@greenwood1967elastic] can be significantly larger than the one predicted by Hertzian theory in case of light indentations, therefore the separation between the micro and macro scales cannot be assumed. The apparent contact area is thus not easy to determine. In theory, because of the infinite support of the Gaussian distribution, the contact has a non-zero probability to occur at any distance $r < R$ from the contact centre. In case of deterministic multi-asperity model, the contact spots appear at random locations outside the Hertzian contact zone. Since the apparent contact area $A_0$ is unknown, the mean contact pressure $\bar p = N/A_0$ is also unknown as well as the contact area fraction $A/A_0$. However, the latter can be evaluated locally with our numerical model by using the stored coarse grained data. The striking difference between the theoretical prediction  and the numerical results , comes from the fact that the former assumes non-interacting asperities and also assumed the scale separation and the validity of Hertz’s theory at macro-scale, which visibly does not hold.

The coarse grained and smoothed with bi-cubic Bézier interpolation results for the spatial area distribution are presented in Fig. \[fig:map\_area\] and compared with Hertzian contact radius for the same force. These results are obtained by averaging the coarse grained results over 2000 simulations with $R=2.5$ mm, $\sigma=10$ $\mu$m, $L/\lambda_l=512$, $L/\lambda_s=24$, $H=0.8$ at different indentations depth ranging from $Z_0/\sigma=4$ to $Z_0/\sigma\approx1.4$, the latter corresponds to the contact force $N\approx116.6$ (N). At very light loads the expansion of the contact zone is remarkably bigger than the Hertzian prediction. This discrepancy decreases under increasing load, but still remains big in the studied loading interval.

![\[fig:map\_area\]Evolution of the contact area fraction under the spherical indentation: the brighter colours correspond to bigger contact area fraction: sub-figures (a-f) correspond to sphere’s tip position $Z_0/\sigma \approx \{4,\;3.57,\;3.14,\;2.71,\;2.29,\;1.86\}$, respectively. The black square is of size $1\times1$ mm. The white dashed line denotes the outline of the Hertzian contact area’s prediction for the mean corresponding contact force $N$.](map_area){width="100.00000%"}

To compute the radial distribution of the contact area fraction and contact pressure at different loads, we use the coarse grained data : $A'_{ij} = A'(x_i,y_j)$, $p_{ij} = p(x_i,y_j)$, where $i,j \in [0,N-1]$ and $x_i = (i+0.5)dx$, $y_j = (j+0.5)dx$. In our simulations, the following coarse-graining was used $N=50$ and $dx=20$ $\mu$m. To find the average radial data, we first interpolated $A'_{ij}$ and $p_{ij}$ in between grid points using bi-cubic Bézier interpolation to obtain $\tilde A'(x,y)$, $\tilde p(x,y)$, then we averaged over concentric circles centred at the indenter’s tip $X_0,Y_0$ to obtain mean radial values $A(r),p'(r)$, which assumes that in statistical sense the both quantities should be axisymmetric: $$A'(r) = \frac{1}{2\pi}\int\limits_0^{2\pi} \tilde A'(X_0+r\cos(\phi),Y_0+r\sin(\phi)) \,d\phi,$$ $$p(r) = \frac{1}{2\pi}\int\limits_0^{2\pi} \tilde p(X_0+r\cos(\phi),Y_0+r\sin(\phi)) \,d\phi.$$ The resulting radial pressure distributions at different load levels are presented in Fig. \[fig:radial\_pres\](a,b) in normal and semi-log scales, respectively. An expected resemblance of our results with those presented in Fig.5(a) from Greenwood and Tripp’s analysis [@greenwood1967elastic] can be noticed. The pressure distribution resembles Gaussian one, but it decays faster than the Gaussian after approximately two times standard deviation from the centre. The extension of the contact zone beyond the Hertzian contact is even more evident in the figures of radial distributions. The contact area fraction which follows approximately the same decay is not shown. It could be noticed that for selected contact characteristics the contact area fraction does not go beyond 1% even at the highest considered load $N\approx116.6$ (N). Therefore the use of the asperity based model can be justified.

![\[fig:radial\_pres\]Contact pressure distribution in rough contact for different loads compared with Hertzian solution (dashed curves of the same colour) for the same contact force (a) normal scale, (b) semi-log scale. Black solid lines correspond to Gaussian fit.](Pres "fig:"){width="100.00000%"}\
(a)\
![\[fig:radial\_pres\]Contact pressure distribution in rough contact for different loads compared with Hertzian solution (dashed curves of the same colour) for the same contact force (a) normal scale, (b) semi-log scale. Black solid lines correspond to Gaussian fit.](Pres_semilog "fig:"){width="100.00000%"}\
(b)

Probability distributions
-------------------------

The deformation regime of individual asperities can be determined by the mean contact pressure at every asperity $\bar p_k = F_k/(\pi a_k^2)$. Therefore it could be useful to construct the probability density function of mean contact pressures $P(\bar p)$, which are experienced at individual asperities. To do so, we carried out simulations for different indentation depths $Z_0/\sigma = \{3,\, 2.5,\, 2\}$ with $R=2.5$ mm, $\sigma=10$ $\mu$m, $L/\lambda_l=512$, $L/\lambda_s=24$, $H=0.8$; for every distribution, the data of 2000 simulations were used, which results in $\approx 320\,000$, $157\,000$, $39\,000$ data points for $Z_0/\sigma = \{2,\, 2.5,\, 3\}$ indentation depths, respectively. The results for $Z_0/\sigma = 2$ are shown in Fig. \[fig:pres\_probability\], the data look the same for other indentation depths. The numerical data fit very accurately the Rayleigh distribution: $$P_R(\bar p) = \frac{\bar p}{\Delta^2}\exp\left(-\frac{\bar p^2}{2\Delta^2}\right),
\label{eq:R}$$ where $\Delta$ is the scaling parameter of pressure units. Clearly, the pressures at asperities are too huge to retain the deformation within the elastic regime. Normally, the maximal contact pressure in the interface is determined by material hardness $H$, which for elasto-plastic materials can be approximately assumed to be $H\approx3\sigma_Y$, where $\sigma_Y$ is the plastic yield stress [@tabor2000hardness; @johnson1985b; @hill1989theoretical; @mesarovic1999spherical; @kogut2002elastic]. The yield stress for gold ranges between $100$-$200$ MPa and the Brinnel hardness ranges in $190-250$ MPa. However, at the scale of asperities the indentation size effect can be pronounced. For spherical indentation, it can be taken into account in terms of hardness change with the indenter’s radius $R$ [@swadener2002correlation]: $$H = H_0\sqrt{1+\frac{R^*}{R}},$$ where $R^* = \bar r/(\rho_s b)$, and $\bar r$ is the Nye factor[^10], $\rho_s$ is the density of statistically stored dislocations, and $b$ is the Burger’s vector. The Burger’s vector in FCC gold with lattice constant $a=0.408$ nm, is given by $b=a/\sqrt2 \approx 0.288$ nm, the density of statistically stored dislocations in strongly deformed region under the indentation could be estimated to be $\rho_s \approx 10^{16}$ m$^{-2}$, taking $\bar r = 2$, we obtain $R^* \approx 0.69$ $\mu$m, which is in good agreement with the data presented in [@kim2018indentation]. Assuming that the mean asperity radius is $\bar r_g \approx 0.737/\sqrt{\bar m_4}$ and assuming macroscopic hardness $H_0=250$ MPa, we obtain the micro-asperity hardness to be $$H = 250 \sqrt{1 + 0.69/\bar r_g}  $$ with $\bar r_g$ measured in micrometers. For the given example with $L/\lambda_l=24$, $L/\lambda_s=512$ and $H=0.8$ and $L=1$ mm, we obtain $\bar r_g \approx 0.18$ $\mu$m, raising the micro-hardness to $H\approx 1.2$ GPa, which is still very small compared to the pressure range going up to 0.5 TPa. Definitely, almost all contacting asperities would be deformed plastically, therefore including elasto-plastic asperity behaviour would be needed for a predictive modelling.

![\[fig:pres\_probability\]Mean asperity pressure probability $\bar p$ evaluated over 2000 simulations for $Z_0/\sigma = 2$ shown in normal and log-scale (see inset). Black lines correspond to fitted Rayleigh distribution .](Pressure_probability_new){width="100.00000%"}

The full pressure distribution can be obtained from this distribution. At every asperity, the Hertzian pressure distribution results in a linear probability density $P_a(p)$: $$P_a(p) = 
  \left\{ \begin{array}{ll}
          \displaystyle 2p/p_0^2 ,&\mbox{ if } p \le p_0\\
          0,&\mbox{ otherwise},
            \end{array}\right.
            \label{eq:asp_pdf}$$ where $p_0$ is the maximal pressure at this asperity, which in turn is given by $p_0 = 1.5\bar p$. The resulting full pressure distribution is shown in Fig. \[fig:full\_distrib\]. Evidently, in the limit of zero pressure the global contact pressure distribution should be linear as it will present a sum of linear contribution from few individual asperities . However, because of the abrupt decay of these probability densities at given $p=p_0$, the global probability density has a slower growth than linear one in the Rayleigh distribution, the decay is also slower than Gaussian. Our simulations suggest, that the global probability distribution is given by the following function: $$P_{gl}(p) = \gamma p^\alpha \exp\left[-(p/\delta)^\beta\right],
\label{eq:my_distr}$$ where we identified $\alpha\approx2/3$ and $\beta\approx3/2$. The fact that $\alpha \approx 2/3$ is consistent with recent numerical findings [@wang2017gauging] where the author identified $\alpha\approx0.7$, which at the same contradicts the Persson’s model which predicts a linear slope for relatively “small” pressures. The sub-Gaussian decay for high pressures can probably be explained by the fact that we consider contact between a sphere and flat rough surface and not between nominally flat surfaces as it is usually the case for the rough contact analyses.

![\[fig:full\_distrib\]Probability density of the total pressure distribution evaluated over 2000 simulations at $Z_0/\sigma=2$. Black lines correspond to Eq.  with $\alpha=2/3$ and $\beta=3/2$.](Total_pres_distr){width="100.00000%"}

The distribution of individual contact spot sizes associated with contacting asperities is of big importance for adhesive wear [@bowden_tabor; @frerot2018mechanistic] as well as for electric and thermal conductivity [@greenwood1966constriction; @slade2017electrical]. The distribution of areas evaluated for same parameters and for indentation depth $Z_0/\sigma=2$ is plotted in Fig. \[fig:area\_distr\] and follows exponential decay $P(A) \sim \exp(-A/a)$. Our results share the same features with experimental observations [@dieterich1996imaging] and numerical results [@hyun2007elastic; @muser2017meeting; @frerot2018mechanistic], namely a “plateau” at small areas in log-log scale and a faster decay at larger areas. However, we do not observe a power law decay in a considerable interval of areas, as was observed in all aforementioned studies. The decay is purely exponential. It can be justified by the fact that individual contact clusters associated with separate asperities never merge in our model, therefore the number of bigger clusters is lower that what is observed in the literature for stronger loads and areas, which considerably overpass $\approx0.5$ % of contact area fraction obtained at these loading conditions with our model.

![\[fig:area\_distr\]Probability density of contact areas associated with individual asperities (normal scale and log-log in the inset), the distribution evaluated over 2000 simulations at $Z_0/\sigma=2$. Black lines correspond to exponential decay $P(A) \sim \exp(-A/a)$](New_area_probability){width="100.00000%"}

\[sec:4\]Discussion
===================

The deterministic multi-asperity model inspired from [@afferrante2012interacting] enables to solve precisely the mechanical normal contact of two spheres in presence of roughness. The semi-analytical solution [@greenwood1967elastic] of this problem was derived by Greenwood and Tripp in late 60’s right after the pioneering work of Greenwood and Williamson on mechanical contact of rough elastic half-spaces [@greenwood1966prcl]. In their solution, the elastic interaction between asperities is taken into account in a statistical sense via an iterative numerical scheme. Recently, Pastewka and Robbins revisited this problem using deterministic full-scale simulations of elastic contact between a sphere and a rough half-space [@pastewka2016contact]. An important conclusion of the both studies is that there’s a linear regime of evolution of the contact area with the applied force. This regime, rather obvious from the physical point of view, extends over many orders of magnitude of the normal load. The added value of Pastewka and Robbins was identification of purely Hertzian contact regime at high loads, and an important extension to adhesive case, which is relevant for nano-metric contact systems.

In our deterministic multi-asperity model with elastic interacting we can properly capture two contact regimes occurring at very light and at moderate loads. In the first regime [@greenwood1970contact], the contact occurs only at a single asperity and thus is accurately described by Hertzian theory of contact. The second regime is determined by the contact of multiple asperities, and results in linear evolution of the contact area as discussed in [@greenwood1967elastic; @pastewka2016contact]. Our model allows to accurately evaluate the transition between the two regimes, which can be easily predicted analytically. The transition load depends on the second and fourth spectral moments of the roughness and does not depend on the indenter’s radius at least as long as it remains much bigger than asperity curvature. Contrary to the contribution of Pastewka and Robbins [@pastewka2016contact], we consider only light indentation depths of order of few surface’s RMS, therefore we operate far from the purely Hertzian regime. That is why, the apparent contact radius is significantly bigger than the one predicted by the Hertzian contact, and the mean pressure, needed to estimate the nominal contact area cannot be properly evaluated. These results once again highlight the lack of scale separation between a macro-scale shape and surface roughness in non-conformal contacts, i.e. rigorously the effect of surface roughness cannot be simply included as an interfacial constitutive law in macroscopic simulations made with smooth surfaces.

The deterministic multi-asperity model can be easily extended to account for non-linear deformation of contacting asperities (viscous and plastic models can be used) while keeping the long-range interaction elastic. Adhesive contact of individual asperities also can be taken into account using either DMT [@DMT], JKR [@JKR] or Maugis-Dugdale model [@maugis1992adhesion; @johnson1997adhesion]. The considered deterministic model does not inherit the main limitation of multi-asperity models, namely the absence of long-range elastic interactions between asperities. The coalescence of contacting asperities [@greenwood2007w; @afferrante2012interacting] in the considered limit of light loads apparently does not present a problem for the considered case, as for realistic root mean squared roughness, the contact area remains below $1$ % in most considered cases.

Statistics of contact pressures clearly confirms the known fact that all contacting asperities operate in severe plastic regime. As was demonstrated here, including the size indentation effect does not change this conclusion. At light loads, the probability distribution of contact areas is found to follow exponential decay. In our model, the total pressure distribution for light pressures follows a power law with exponent $2/3$ and at higher loads decays as $\exp\left(-(p/p_0)^{3/2}\right)$, which is different from what Persson’s model predicts for two nominally flat surfaces: linear growth of probability at light pressures interval, and a Gaussian decay for higher pressures. Nevertheless, the demonstrated sub-linear growth of the pressure probability of relatively small pressure is in good agreement with recent results obtained with full deterministic numerical simulations [@wang2017gauging].

[**Acknowledgment**]{}\

Helpful comments and remarks of Andrei Shvarts on the paper are greatly appreciated.

[**Appendix**]{}

Spectral density and asperity density\[app\]
============================================

Consider a surface spectrum $\Phi$ with a plateau which corresponds to an isotropic and periodic rough surface with a period $L$. Following the construction technique [@hu1992ijmtm] for all possible directions determined by wavenumbers $k_x,k_y$ in $OX$ and $OY$ directions, respectively, in average for $K=\sqrt{k_x^2+k_y^2}$, the following equation holds: $$\langle \Phi(K) \rangle =\left\{ \begin{array}{ll}
          \Phi_0 ,&\mbox{ if } \xi \le K/k_l \le 1\\
          \Phi_0 (K/k_l)^{-2(1+H)},&\mbox{ if } 1 < K/k_l \le \zeta\\
          0,&\mbox{ otherwise},
            \end{array}\right.$$ where $k_l = 2\pi/\lambda_l$, $k_s=2\pi/\lambda_s$, $\xi = 2\pi/(Lk_l)$, $\zeta = \lambda_l/\lambda_s = k_s/k_l$, for the following inequalities $\lambda_s \le \lambda_l \le L$, where from it follows $\xi \le 1$, $\zeta \ge 1$. A spectral moment $m_p$ is given by the following convolution $$m_{p} = \int\limits_{\frac{2\pi}{L}}^{k_s} K^{p+1} \Phi(K) dK \int\limits_{0}^{2\pi}\cos^{p}(\varphi)d\varphi  =
 \Phi_0T(p)\left[\left.\frac{K^{p+2}}{p+2}\right|_{\frac{2\pi}{L}}^{k_l} + \left.\frac{K^{p-2H}k_l^{2(1+H)}}{(p-2H)}\right|_{k_l}^{k_s}\right]$$ $$m_{p} = \Phi_0 k_l^{p+2} T(p)\left[\frac{1-\xi^{p+2} }{p+2} + \frac{\zeta^{p-2H} - 1 }{p-2H}\right],
 \label{eq:a:mq0}$$ where $T(p)$ $$T(p) = \int\limits_{0}^{2\pi} \cos^p(\varphi)d\varphi,\quad T(0) = 2\pi, \quad T(2) = \pi,\quad T(4) = 3\pi/4.$$ The density of asperities can be found as [@nayak1971tasme] $$D = \frac{\sqrt3}{18\pi}\frac{m_4}{m_2},$$ which by substituting expression for spectral moments  takes the following form [@yastrebov2015ijss]: $$\displaystyle
D = \frac{\sqrt3 k_l^2}{24\pi} \left(\frac{1-\xi^{6}}{3} + \frac{\zeta^{4-2H} - 1 }{2-H}\right) \left/\left(\frac{1 - \xi^{4}}{2} + \frac{\zeta^{2-2H} - 1 }{1-H}\right)\right.,$$ which for a surface without a plateau in the spectral density takes a simpler form $$\displaystyle
D \approx \frac{\sqrt3 k_l^2}{24\pi} \frac{1-H}{2-H}\frac{\zeta^{4-2H}}{\zeta^{2-2H} - 1 }.
\label{eq:density}$$ Note that in the numerator we readily approximated $(\zeta^{4-2H}-1)$ by $\zeta^{4-2H}$; however, this approximation cannot be made in the denominator if $H$ is reasonably high. In Fig. \[fig:density\] we compare this analytical expression  with numerically computed asperity densities for surfaces with different Hurst exponents $H \in [0.1,0.9]$ and different magnifications $\zeta = k_s/k_l \in [3.125, 125]$, to keep the surfaces representative we fixed the longest wavelength to $L/\lambda_l = 8$ and considered surfaces without a plateau. Note that to ensure accuracy of asperity detection, for $L/\lambda \le 500$ ($\zeta \le 62.5$) we used surfaces with discretization $N = 4096$, and for $500 < L/\lambda_s \le 1000$ ($62.5  < \zeta \le  125$), discretization $N = 8192$ was used. This choice was made because a considerable underestimation of asperity densities was observed for insufficient surface discretization, e.g. the relative error for $N=4096$, $L/\lambda_l = 8$, $L/\lambda_s = 1000$, reached $\approx 8$% for $H=0.9$ and raised to $\approx 15\%$ for $H=0.1$.

![\[fig:density\]Comparison of Eq.  with numerically evaluated asperity density of generated rough surfaces with $L/\lambda_l = 8$. Cross points of the analytical solution correspond to the same parameters used in numerical estimation.](Density_asperities){width="100.00000%"}

The average distance between closest asperities for a surface without plateau is given by: $$\langle d\rangle = 1/\sqrt{D} \approx \sqrt{\frac{24\pi}{\sqrt3 k_l^2} \frac{(2-H)}{(1-H)}\frac{(\zeta^{2-2H} - 1)}{\zeta^{4-2H}}} = \lambda_l \sqrt{\frac{6 }{\sqrt3  \pi } \frac{(2-H)}{(1-H)}\frac{(\zeta^{2-2H} - 1)}{\zeta^{4-2H}}}.$$ As close asperities have a comparable height, the critical contact radius for asperities, at which the associated contact zones coalesce, may be estimated as $a' = \langle d\rangle/2$. Moreover, at average contact radius $a_c = \langle d\rangle/4$ the underlying surface would change curvature and thus the Hertzian assumption cannot be valid anymore. Therefore we shall assume that that the critical contact radius for asperities, beyond which the solution is not valid, is given by $$a_c = \lambda_l \sqrt{\frac{\sqrt3 }{8  \pi } \frac{(2-H)}{(1-H)}\frac{(\zeta^{2-2H} - 1)}{\zeta^{4-2H}}}.$$

[10]{}

L. Afferrante, F. Bottiglione, C. Putignano, B. Persson, and G. Carbone. Elastic contact mechanics of randomly rough surfaces: An assessment of advanced asperity models and persson’s theory. , 66(2):75, 2018.

L. Afferrante, G. Carbone, and G. Demelio. Interacting and coalescing hertzian asperities: A new multiasperity contact model. , 278:28–33, 2012.

L. Afferrante, G. Carbone, and G. Demelio. Interacting and coalescing hertzian asperities: a new multiasperity contact model. , 278:28–33, 2012.

R. Aghababaei, D. H. Warner, and J.-F. Molinari. Critical length scale controls adhesive wear mechanisms. , 7, 2016.

S. Akarapu, T. Sharp, and M. O. Robbins. Stiffness of contacts between rough surfaces. , 106(20):204301, 2011.

G. Anciaux and J. F. Molinari. Contact mechanics at the nanoscale, a 3d multiscale approach. , 79(9):1041–1067, 2009.

J. Archard and W. Hirst. The wear of metals under unlubricated conditions. , 236(1206):397–410, 1956.

J. Barber. Bounds on the electrical resistance between contacting elastic rough bodies. In [*Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*]{}, volume 459, pages 53–66. The Royal Society, 2003.

S. Basseville, E. H[é]{}ripr[é]{}, and G. Cailletaud. Numerical simulation of the third body in fretting problems. , 270(11-12):876–887, 2011.

A. Bemporad and M. Paggi. Optimization algorithms for the solution of the frictionless normal contact between rough surfaces. , 69:94–105, 2015.

O. Ben-David, G. Cohen, and J. Fineberg. The dynamics of the onset of frictional slip. , 330(6001):211–214, 2010.

F. P. Bowden and D. Tabor. . Oxford University Press, 2001.

A. W. Bush, R. D. Gibson, and T. R. Thomas. The elastic contact of a rough surface. , 35(1):87–111, 1975.

C. Campan[á]{} and M. H. M[ü]{}ser. Practical green’s function approach to the simulation of elastic semi-infinite solids. , 74(7):075420, 2006.

G. Carbone and F. Bottiglione. Asperity contact theories: Do they predict linearity between contact area and load? , 56:2555–2572, 2008.

M. Ciavarella. On the significance of asperity models predictions of rough contact with respect to recent alternative theories. , 139(2):021402, 2017.

W. B. Dapp, A. L[ü]{}cke, B. N. Persson, and M. H. M[ü]{}ser. Self-affine elastic contacts: percolation and leakage. , 108(24):244301, 2012.

W. B. Dapp and M. H. M[ü]{}ser. Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions. , 109(4):44001, 2015.

W. B. Dapp and M. H. M[ü]{}ser. Fluid leakage near the percolation threshold. , 6, 2016.

B. V. Derjaguin, V. M. Muller, and Y. P. Toporov. Effect of contact deformations on the adhesion of particles. , 53(2):314–326, 1975.

T. Dick and G. Cailletaud. Fretting modelling with a crystal plasticity model of ti6al4v. , 38(1):113–125, 2006.

J. H. Dieterich and B. D. Kilgore. Imaging surface contacts: power law contact distributions and contact stresses in quartz, calcite, glass and acrylic plastic. , 256(1-4):219–239, 1996.

S. Fouvry, P. Kapsa, H. Zahouani, and L. Vincent. Wear analysis in fretting of hard coatings through a dissipated energy concept. , 203:393–403, 1997.

L. Fr[é]{}rot, R. Aghababaei, and J.-F. Molinari. A mechanistic understanding of the wear coefficient: From single to multiple asperities contact. , 114:172–184, 2018.

K. Fuller and D. Tabor. The effect of surface roughness on the adhesion of elastic solids. , 345(1642):327–342, 1975.

I. Goryacheva. Mechanics of discrete contact. , 39(5):381–386, 2006.

I. G. Goryacheva. , volume 61. Springer Science & Business Media, 2013.

J. Greenwood and J. Tripp. The contact of two nominally flat rough surfaces. , 185(1):625–633, 1970.

J. A. Greenwood. Constriction resistance and the real area of contact. , 17(12):1621, 1966.

J. A. Greenwood. Formulas for moderately elliptical hertzian contacts. , 107:501–504, 1985.

J. A. Greenwood. A simplified elliptic model of rough surface contact. , 261:191–200, 2006.

J. A. Greenwood. A note on [N]{}ayak’s third paper. , 262:225–227, 2007.

J. A. Greenwood and J. H. Tripp. The elastic contact of rough spheres. , 34(1):153–159, 1967.

J. A. Greenwood and J. B. P. Williamson. Contact of nominally flat surfaces. , 295:300–319, 1966.

R. Hill, B. Stor[å]{}kers, and A. Zdunek. A theoretical study of the brinell hardness test. , 423(1865):301–330, 1989.

D. A. Hills. Mechanics of fretting fatigue. , 175(1-2):107–113, 1994.

K. Hokkirigawa and K. Kato. An experimental and theoretical investigation of ploughing, cutting and wedge formation during abrasive wear. , 21(1):51 – 57, 1988.

Y. Z. Hu and K. Tonder. Simulation of [3-D]{} random rough surface by [2-D]{} digital filter and fourier analysis. , 32:83–90, 1992.

G. Huber, S. N. Gorb, N. Hosoda, R. Spolenak, and E. Arzt. Influence of surface roughness on gecko adhesion. , 3(4):607–610, 2007.

S. Hyun, L. Pei, J. F. Molinari, and M. O. Robbins. Finite-element analysis of contact between elastic self-affine surfaces. , 70(2):026117, 2004.

S. Hyun and M. O. Robbins. Elastic contact between rough surfaces: Effect of roughness at large and small wavelengths. , 40(10-12):1413–1422, 2007.

K. Johnson and J. Greenwood. An adhesion map for the contact of elastic spheres. , 192(2):326–333, 1997.

K. L. Johnson. . Cambridge University Press, 1985.

K. L. Johnson, K. Kendall, and A. Roberts. Surface energy and the contact of elastic solids. , 324(1558):301–313, 1971.

Y.-C. Kim, E.-J. Gwak, S.-m. Ahn, N.-R. Kang, H. N. Han, J.-i. Jang, and J.-Y. Kim. Indentation size effect for spherical nanoindentation on nanoporous gold. , 143:10–14, 2018.

L. Kogut and I. Etsion. Elastic-plastic contact analysis of a sphere and a rigid flat. , 69(5):657–662, 2002.

K. Kubiak, T. Liskiewicz, and T. Mathia. Surface morphology in engineering applications: Influence of roughness on sliding and wear in dry fretting. , 44(11):1427 – 1432, 2011.

C. V. Madhusudana. . Springer, 1996.

D. Maugis. Adhesion of spheres: the jkr-dmt transition using a dugdale model. , 150(1):243–269, 1992.

P. Meakin. , volume 5. Cambridge university press, 1998.

S. D. Mesarovic and N. A. Fleck. Spherical indentation of elastic–plastic solids. In [*Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*]{}, volume 455, pages 2707–2728. The Royal Society, 1999.

Y. Mo, K. T. Turner, and I. Szlufarska. Friction laws at the nanoscale. , 457(7233):1116, 2009.

M. H. M[ü]{}ser, W. B. Dapp, R. Bugnicourt, P. Sainsot, N. Lesaffre, T. A. Lubrecht, B. N. Persson, K. Harris, A. Bennett, K. Schulze, et al. Meeting the contact-mechanics challenge. , 65(4):118, 2017.

P. R. Nayak. Random process model of rough surfaces. , 93:398–407, 1971.

M. Paggi and M. Ciavarella. The coefficient of proportionality $\kappa$ between real contact area and load, with new asperity models. , 268(7-8):1020–1029, 2010.

L. Pastewka, S. Moser, P. Gumbsch, and M. Moseler. Anisotropic mechanical amorphization drives wear in diamond. , 10(1):34, 2011.

L. Pastewka, N. Prodanov, B. Lorenz, M. H. M[ü]{}ser, M. O. Robbins, and B. N. J. Persson. Finite-size scaling in the interfacial stiffness of rough elastic contacts. , 87:062809, Jun 2013.

L. Pastewka and M. O. Robbins. Contact between rough surfaces and a criterion for macroscopic adhesion. , page 201320846, 2014.

L. Pastewka and M. O. Robbins. Contact area of rough spheres: large scale simulations and simple scaling laws. , 108(22):221601, 2016.

N. Patir and H. Cheng. An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication. , 100(1):12–17, 1978.

L. Pei, S. Hyun, J. Molinari, and M. O. Robbins. Finite element modeling of elasto-plastic contact between rough surfaces. , 53(11):2385–2409, 2005.

L. Pei, S. Hyun, J. F. Molinari, and M. O. Robbins. Finite element modeling of elasto-plastic contact between rough surfaces. , 53:2385–2409, 2005.

B. N. J. Persson. Theory of rubber friction and contact mechanics. , 115:3840–3861, 2001.

I. Polonsky and L. Keer. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques. , 231(2):206–219, 1999.

H. Proudhon, S. Fouvry, and J.-Y. Buffière. A fretting crack initiation prediction taking into account the surface roughness and the crack nucleation process volume. , 27(5):569 – 579, 2005.

C. Putignano, L. Afferrante, G. Carbone, and G. Demelio. The influence of the statistical properties of self-affine surfaces in elastic contacts: A numerical investigation. , 60:973–982, 2012.

C. Putignano, L. Afferrante, G. Carbone, and G. Demelio. A new efficient numerical method for contact mechanics of rough surfaces. , 49(2):338–343, 2012.

A. Ramesh and S. N. Melkote. Modeling of white layer formation under thermally dominant conditions in orthogonal machining of hardened aisi 52100 steel. , 48(3-4):402–414, 2008.

V. Rey, S. Krumscheid, and F. Nobile. Quantifying uncertainties in contact mechanics of rough surfaces using the multilevel monte carlo method. , 2017. available online on hal.archives-ouvertes.fr/hal-01596319/.

J. R. Rice. Heating and weakening of faults during earthquake slip. , 111(B5), 2006.

A. Ruina. Slip instability and state variable friction laws. , 88(B12):10359–10370, 1983.

T. Sasada, M. Oike, and N. Emori. The effect of abrasive grain size on the transition between abrasive and adhesive wear. , 97(3):291 – 302, 1984.

R. Sayles and P. Macpherson. Influence of wear debris on rolling contact fatigue. In [*Rolling contact fatigue testing of bearing steels*]{}. ASTM International, 1982.

M. Scherge, S. Gorb, and S. N. Gorb. . Springer Science & Business Media, 2001.

A. G. Shvarts and V. A. Yastrebov. Trapped fluid in contact interface. , 2018.

P. G. Slade. . CRC press, 2017.

H. M. Stanley and T. Kato. An [FFT]{}-based method for rough surface contact. , 119:481–485, 1997.

J. Swadener, E. George, and G. Pharr. The correlation of the indentation size effect measured with indenters of various shapes. , 50(4):681–694, 2002.

D. Tabor. . Oxford university press, 2000.

A. Vakis, V. Yastrebov, J. Scheibert, L. Nicola, D. Dini, C. Minfray, A. Almqvist, M. Paggi, S. Lee, G. Limbert, J. Molinari, G. Anciaux, R. Aghababaei, S. E. Restrepo, A. Papangelo, A. Cammarata, P. Nicolini, C. Putignano, G. Carbone, S. Stupkiewicz, J. Lengiewicz, G. Costagliola, F. Bosia, R. Guarino, N. Pugno, M. M[ü]{}ser, and M. Ciavarella. Modeling and simulation in tribology across scales: An overview. , 125:169 – 199, 2018.

A. Wang and M. H. M[ü]{}ser. Gauging persson theory on adhesion. , 65(3):103, 2017.

P. Wriggers and G. Zavarise. Application of augmented lagrangian techniques for non-linear constitutive laws in contact interfaces. , 9(10):815–824, 1993.

V. Yastrebov, G. Anciaux, and J.-F. Molinari. The contact of elastic regular wavy surfaces revisited. , 56(1):171–183, 2014.

V. A. Yastrebov, G. Anciaux, and J.-F. Molinari. Contact between representative rough surfaces. , 86(3):035601(R), 2012.

V. A. Yastrebov, G. Anciaux, and J.-F. Molinari. From infinitesimal to full contact between rough surfaces: evolution of the contact area. , 52:83–102, 2015.

V. A. Yastrebov, G. Anciaux, and J.-F. Molinari. The role of the roughness spectral breadth in elastic contact of rough surfaces. , 107:469–493, 2017.

V. A. Yastrebov, G. Cailletaud, H. Proudhon, F. S. M. Mballa, S. No[ë]{}l, P. Test[é]{}, and F. Houz[é]{}. Three-level multi-scale modeling of electrical contacts sensitivity study and experimental validation. In [*Electrical Contacts (Holm), 2015 IEEE 61st Holm Conference on*]{}, pages 414–422. IEEE, 2015.

V. A. Yastrebov, J. Durand, H. Proudhon, and G. Cailletaud. Rough surface contact analysis by means of the finite element method and of a new reduced model. , 339(7-8):473–490, 2011.

V. A. Yastrebov, J. Durand, H. Proudhon, and G. Cailletaud. Rough surface contact analysis by means of the finite element method and of a new reduced model. , 339:473–490, 2011.

X. Yin and K. Komvopoulos. An adhesive wear model of fractal surfaces in normal contact. , 47(7-8):912–921, 2010.

[^1]: We recall that in the vicinity of the tip, a sphere and a paraboloid of revolution of the same curvature are indistinguishable.

[^2]: It can be assumed that $z$-curvature of an in-plane line connecting two neighbouring asperities, changes at quarter distance from each summit.

[^3]: l.u. is length units.

[^4]: In the original paper [@greenwood2006w], the probability density should be divided by a factor of two to fulfil the normalization.

[^5]: By separation we understand the distance between the sphere’s tip and the mean line of the rough surface.

[^6]: If one assumes that the elastic asperities readily plastifies, the linearity can be preserved even at the first contact, since the contact pressure readily saturates at the hardness level and thus the contact area grows proportionally to the contact force.

[^7]: A more accurate estimation only for highest asperities could be obtained using Eq. (11) from [@greenwood2006w] and by correcting it with a factor $1/2$.

[^8]: Lower index “PR” was used to indicate the authors (Pastewka & Robbins) of the mentioned paper [@pastewka2016contact], we show here the combination of Eq.(5) and (6) from this paper.

[^9]: This scaling factor is implicitly contained in authors’ [@pastewka2016contact] estimation of the mean contact pressure $N/A_{\mbox{\tiny Hertz}}$, which disappears in the transition to the fully Hertzian regime at high loads, but had to be preserved at light loads, which was not the case.

[^10]: The Nye factor is a dimensionless factor determining how much greater dislocations should be created to accommodate plastic deformation compared to the number of geometrically necessary dislocations
