Vladislav A Yastrebov

The Elastic Contact of Rough Spheres Investigated Using a Deterministic Multi-Asperity Model

Keywords: roughness, indentation, contact, deterministic multi-asperity model, contact area

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.

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, contactrelated mechanisms (friction, adhesion, and wear) are still not fully understood and thus are among the most cutting-edge research topics in the mechanical community [START_REF] Mo | Friction laws at the nanoscale[END_REF][START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF][START_REF] Pastewka | Anisotropic mechanical amorphization drives wear in diamond[END_REF][START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF].

Numerous models of contact-related mechanisms exist at structural scale [START_REF] Vakis | Modeling and simulation in tribology across scales: An overview[END_REF]. They serve to model (1) interfacial normal and tangential stiffness [START_REF] Wriggers | Application of augmented lagrangian techniques for non-linear constitutive laws in contact interfaces[END_REF][START_REF] Akarapu | Stiffness of contacts between rough surfaces[END_REF], [START_REF] Afferrante | Interacting and coalescing hertzian asperities: A new multiasperity contact model[END_REF] frictional resistance [START_REF] Ruina | Slip instability and state variable friction laws[END_REF], (3) material removal on rubbing surfaces (wear) [START_REF] Archard | The wear of metals under unlubricated conditions[END_REF][START_REF] Yin | An adhesive wear model of fractal surfaces in normal contact[END_REF], [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF] heat transfer between contacting solids [START_REF] Madhusudana | Thermal contact conductance[END_REF], [START_REF] Akarapu | Stiffness of contacts between rough surfaces[END_REF] contact electric resistance [START_REF] Slade | Electrical contacts: principles and applications[END_REF][START_REF] Yastrebov | Threelevel multi-scale modeling of electrical contacts sensitivity study and experimental validation[END_REF], [START_REF] Anciaux | Contact mechanics at the nanoscale, a 3d multiscale approach[END_REF] adhesion [START_REF] Fuller | The effect of surface roughness on the adhesion of elastic solids[END_REF][START_REF] Pastewka | Contact between rough surfaces and a criterion for macroscopic adhesion[END_REF], [START_REF] Archard | The wear of metals under unlubricated conditions[END_REF] interfacial fluid flow [START_REF] Dapp | Self-affine elastic contacts: percolation and leakage[END_REF][START_REF] Dapp | Fluid leakage near the percolation threshold[END_REF][START_REF] Shvarts | Trapped fluid in contact interface[END_REF], [START_REF] Barber | Bounds on the electrical resistance between contacting elastic rough bodies[END_REF] microstructural changes in near-contact material layers [START_REF] Ramesh | Modeling of white layer formation under thermally dominant conditions in orthogonal machining of hardened aisi 52100 steel[END_REF], [START_REF] Basseville | Numerical simulation of the third body in fretting problems[END_REF] fretting wear life-cycle [START_REF] Hills | Mechanics of fretting fatigue[END_REF][START_REF] Fouvry | Wear analysis in fretting of hard coatings through a dissipated energy concept[END_REF][START_REF] Dick | Fretting modelling with a crystal plasticity model of ti6al4v[END_REF], [START_REF] Bemporad | Optimization algorithms for the solution of the frictionless normal contact between rough surfaces[END_REF] debris generation [START_REF] Sayles | Influence of wear debris on rolling contact fatigue[END_REF][START_REF] Basseville | Numerical simulation of the third body in fretting problems[END_REF], [START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF] lubrication [START_REF] Patir | An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication[END_REF], 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. [START_REF] Afferrante | Elastic contact mechanics of randomly rough surfaces: An assessment of advanced asperity models and persson's theory[END_REF] 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. [START_REF] Afferrante | Interacting and coalescing hertzian asperities: a new multiasperity contact model[END_REF] In case of adhesive wear, the removed volume depends on the characteristics and number of contact spots forming the true contact area. [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF] 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. [START_REF] Akarapu | Stiffness of contacts between rough surfaces[END_REF] 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. [START_REF] Anciaux | Contact mechanics at the nanoscale, a 3d multiscale approach[END_REF] 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 [START_REF] Scherge | Biological micro-and nanotribology[END_REF][START_REF] Huber | Influence of surface roughness on gecko adhesion[END_REF] is more related to the hierarchical structure of their limbs. [START_REF] Archard | The wear of metals under unlubricated conditions[END_REF] 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. [START_REF] Barber | Bounds on the electrical resistance between contacting elastic rough bodies[END_REF] 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 [START_REF] Bowden | The Friction and Lubrication of Solids[END_REF][START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]. These extreme thermo-mechanical loads can induce phase transformations and changes in material microstructures. [START_REF] Basseville | Numerical simulation of the third body in fretting problems[END_REF] Local wear and initiation of fatigue cracks in fretting conditions [START_REF] Proudhon | A fretting crack initiation prediction taking into account the surface roughness and the crack nucleation process volume[END_REF][START_REF] Kubiak | Surface morphology in engineering applications: Influence of roughness on sliding and wear in dry fretting[END_REF] strongly depend on the local stress state, which is intimately linked with the surface roughness. [START_REF] Bemporad | Optimization algorithms for the solution of the frictionless normal contact between rough surfaces[END_REF] Produced by adhesive wear [START_REF] Sasada | The effect of abrasive grain size on the transition between abrasive and adhesive wear[END_REF], 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 [START_REF] Hokkirigawa | An experimental and theoretical investigation of ploughing, cutting and wedge formation during abrasive wear[END_REF]. The presence of debris (or mobile third body) can significantly affect the true contact area and thus modify the kinetics of debris production. [START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF] 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 [START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF][START_REF] Campaná | Practical green's function approach to the simulation of elastic semi-infinite solids[END_REF][START_REF] Anciaux | Contact mechanics at the nanoscale, a 3d multiscale approach[END_REF][START_REF] Yastrebov | Rough surface contact analysis by means of the finite element method and of a new reduced model[END_REF][START_REF] Putignano | A new efficient numerical method for contact mechanics of rough surfaces[END_REF][START_REF] Afferrante | Interacting and coalescing hertzian asperities: A new multiasperity contact model[END_REF][START_REF] Bemporad | Optimization algorithms for the solution of the frictionless normal contact between rough surfaces[END_REF][START_REF] Dapp | Meeting the contact-mechanics challenge[END_REF] 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 [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF][START_REF] Bush | The elastic contact of a rough surface[END_REF][START_REF] Persson | Theory of rubber friction and contact mechanics[END_REF][START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load?[END_REF] 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 [START_REF] Greenwood | The elastic contact of rough spheres[END_REF][START_REF] Goryacheva | Mechanics of discrete contact[END_REF][START_REF] Goryacheva | Contact mechanics in tribology[END_REF][START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF], we choose to use a deterministic multi-asperity approach, rather similar to ones elaborated in [START_REF] Afferrante | Interacting and coalescing hertzian asperities: a new multiasperity contact model[END_REF][START_REF] Afferrante | Elastic contact mechanics of randomly rough surfaces: An assessment of advanced asperity models and persson's theory[END_REF], 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 [START_REF] Greenwood | The elastic contact of rough spheres[END_REF]. The deterministic multi-asperity model, compared to statistical multi-asperity models [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF][START_REF] Bush | The elastic contact of a rough surface[END_REF][START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load?[END_REF], 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 [START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF][START_REF] Yastrebov | Rough surface contact analysis by means of the finite element method and of a new reduced model[END_REF] or modifications of boundary element methods [START_REF] Stanley | An FFT-based method for rough surface contact[END_REF][START_REF] Polonsky | A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques[END_REF][START_REF] Campaná | Practical green's function approach to the simulation of elastic semi-infinite solids[END_REF][START_REF] Bemporad | Optimization algorithms for the solution of the frictionless normal contact between rough surfaces[END_REF].

The paper is organized as follows. In Section 2 we define the problem, recall the roughness generation method and explain the method to extract asperities' data. In Section 3, the numerical model to solve contact problems is presented. The results are presented in Section 4, 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 5.

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. 1(a)] of curvature radii R 1 and R 2 , the solids are assumed to be elastic, homogeneous and isotropic with elastic properties E 1 , ν 1 and E 2 , ν 2 , respectively, where E and ν are the Young's modulus and Poisson's ratio, respectively. Schematically this problem can be presented as a contact between two spheres1 , see Fig. 1. 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 1 R 2 /(R 1 + R 2 ) and a rough elastic half-space with effective elastic modulus E * [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Barber | Bounds on the electrical resistance between contacting elastic rough bodies[END_REF]:

E * = E 1 E 2 (1 -ν 2 2 )E 1 + (1 -ν 2 1 )E 2 .
In case of smooth surfaces, we recall [START_REF] Johnson | Contact mechanics[END_REF] below the analytical solutions for this problem under the action of a squeezing force N. The contact radius a, the approach between distance points δ, and the pressure distribution p(x, y) are given by the following expressions:

a = 3NR * 4E * 1/3 , δ = a 2 R * = 9N 2 16R * E * 2 1/3 , p(x, y) = 3P 2πa 2 1 - x 2 + y 2 a 2 . ( 1
)

Roughness

To model the surface roughness we use a white-noise filtering technique [START_REF] Hu | Simulation of 3-D random rough surface by 2-D digital filter and fourier analysis[END_REF][START_REF] Yastrebov | The role of the roughness spectral breadth in elastic contact of rough surfaces[END_REF]. 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 λ) of spectral energy Φ ∼ λ 2(1+H) , and also (2) the lower and (3) the upper cut-off wavelengths λ s and λ 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 λ ∈ [λ s , λ l ] can be followed by a plateau for longer wavelengths λ ∈ (λ l , L], where L is the sample length.

R 2 R 1 E 1 ,ν 1 E 2 ,ν 2 z z 1 z 2 R * E * z z * N N N
Based on this definition, we can readily conclude that for loads implying the value of the macroscopic contact radius a macro which is comparable with the shortest roughness wavelength a macro ∼ λ s , the contact will be formed by a few contact spots, whereas for a macro λ 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) λ l determines the dispersion of results [START_REF] Yastrebov | Contact between representative rough surfaces[END_REF][START_REF] Pastewka | Finite-size scaling in the interfacial stiffness of rough elastic contacts[END_REF][START_REF] Yastrebov | From infinitesimal to full contact between rough surfaces: evolution of the contact area[END_REF][START_REF] Rey | Quantifying uncertainties in contact mechanics of rough surfaces using the multilevel monte carlo method[END_REF]: 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 macro λ 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 [START_REF] Afferrante | Interacting and coalescing hertzian asperities: A new multiasperity contact model[END_REF] 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 [START_REF] Yastrebov | The contact of elastic regular wavy surfaces revisited[END_REF][START_REF] Dapp | Contact mechanics of and reynolds flow through saddle points: On the coalescence of contact patches and the leakage rate through near-critical constrictions[END_REF], 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 quarter2 of the average distance between asperities d (see A for details), i.e.

a micro < d 4 = 1 4 √ D ≈ λ l √ 3 8π (2 -H) (1 -H) (ζ 2-2H -1) ζ 4-2H (2)
where D is the density of asperities, and the last approximate equality is valid for surfaces with a rich spectral content, i.e. for ζ = λ l /λ s 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 asperityheight distribution on a semi-infinite support [START_REF] Ciavarella | On the significance of asperity models predictions of rough contact with respect to recent alternative theories[END_REF]. 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 (2). 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 asp predicted by a multi-asperity model over nominal area A 0 is formed by N identical contact spots πa 2 micro , 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 micro ≈ A asp πDA 0 , (3)
where from using (2), a simple limit for the validity of multi-asperity models could be established:

A asp A 0 < π 16 ≈ 19.6%. ( 4
)
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 nonuniform distribution of asperities and contact spots, the real limit can be considerably smaller (see also comparison between deterministic and multi-asperity models [START_REF] Yastrebov | The role of the roughness spectral breadth in elastic contact of rough surfaces[END_REF]).

Asperities

In vicinity of every summit, a smooth surface can be approximated by a quadratic form of local coordinates:

z(x, y) = z 0 + ax 2 + by 2 + 2cxy = z 0 + x y a c c b x y = z 0 + X AX with a = 1 2 ∂ 2 z ∂x 2 , b = 1 2 ∂ 2 z ∂y 2 , c = 1 2 ∂ 2 z ∂x∂y .
To determine the surface principal curvatures and their orientation, we seek for such rotation ϕ of {x, y} coordinates

X = QX , X = x y , X = x y , Q = cos(ϕ) sin(ϕ) -sin(ϕ) cos(ϕ)
that the resulting form is diagonal:

z(x, y) = z 0 + X Q AQX , Q AQ = λ 1 0 0 λ 2 .
The angle ϕ that ensures the diagonality should verify the following equation:

(a -b) cos(ϕ) sin(ϕ) + c(cos 2 (ϕ) -sin 2 (ϕ)) = 0,
which after dividing by (c cos 2 (ϕ)) gives a quadratic equation for tan(ϕ):

tan 2 (ϕ) - a -b c tan(ϕ) -1 = 0,
therefore the two solutions are given by :

ϕ i = arctan          a -b 2c ± a -b 2c 2 + 1          , (5)
which correspond to eigenvalues

λ i = 1 2 a + b ± (a -b) 2 + 4c 2 ,
therefore the principal curvatures and principal radii are given by:

κ i = a + b ± (a -b) 2 + 4c 2 , R i = 1 a + b ± (a -b) 2 + 4c 2 (6)
Rough surface measurements are inevitably of discrete nature, therefore the question of asperity identification is relevant [START_REF] Paggi | The coefficient of proportionality κ between real contact area and load, with new asperity models[END_REF]. Numerically generated rough surfaces [START_REF] Hu | Simulation of 3-D random rough surface by 2-D digital filter and fourier analysis[END_REF][START_REF] Meakin | Fractals, scaling and growth far from equilibrium[END_REF] (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 i j = z(i∆x, j∆y), where ∆x, ∆y are the sampling intervals in x, y directions, respectively. For such a representation, every potential summit z i j is identified if the nearest and next-nearest (diagonal) neighbour points are lower: z i j > z (k + i)∆x, (l + j)∆y , 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 A can be estimated as:

a = 1 2 ∂ 2 z ∂x 2 ≈ z i+1 j + z i-1 j -2z i j 2∆x 2 ; b = 1 2 ∂ 2 z ∂y 2 ≈ z i j+1 + z i j-1 -2z i j 2∆y 2 (7) c = 1 2 ∂ 2 z ∂x∂y ≈ z i+1j+1 + z i-1 j-1 -z i+1j-1 -z i-1 j+1 8∆x∆y . (8)
The principal curvatures and radii then can be found from [START_REF] Anciaux | Contact mechanics at the nanoscale, a 3d multiscale approach[END_REF]. If both curvatures are negative the generalized summit is the real summit. If the principal curvatures have a different sign κ 1 κ 2 < 0 this "summit" is a saddle point masked by discretization. The algebraic κ a and the geometric κ g mean curvatures are defined as follows

κ a = 1 2 (κ 1 + κ 2 ) = 1 2 (a + b) (9)
κ g = √ κ 1 κ 2 = 2 √ ab -c 2 (10)
In Fig. 2 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) (zz) 2 /L = σ/L = 0.025. Obviously, a decrease of the RMS height would result in an increase of asperities' radii. It was shown by Greenwood [START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF], that according to Nayak's theory [START_REF] Nayak | Random process model of rough surfaces[END_REF], the asperities are only "mildly" elliptic, and thus, the geometric mean curvature can be successfully used in multi-asperity contact models [START_REF] Greenwood | Formulas for moderately elliptical hertzian contacts[END_REF] instead of treatment elliptic contacts [START_REF] Bush | The elastic contact of a rough surface[END_REF][START_REF] Johnson | Contact mechanics[END_REF], without any significant loss in accuracy. However, as seen by "naked eye" from Fig. 2, 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 [START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF].

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 [START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF], who obtained the following equation for the joint probability density of the ensemble of asperities' principal curvatures 4 :

P(κ 1 , κ 2 ) = 27 16 √ π κ 1 κ 2 |κ 1 -κ 2 | exp - 3 16 3κ 2 1 + 3κ 2 2 -2κ 1 κ 2 , (11)
where

κ i = κ i / √ m 4 represent normalized curvatures.
This comparison is shown in Fig. 3. Numerically, the statistics is performed over an ensemble of asperities detected on a single realization of a surface with a rich spectrum L/λ l = 4, L/λ s = 2048, L/∆x = 8192 and H = 0.8; in total this surface contains ≈ 580 000 asperities. In accordance with Greenwood's findings, the most probable asperities possess the following ratio of principal curvatures κ 1 /κ 2 ≈ 2.214; the probability to find an asperity close-tospherical is marginal. The overall mean geometrical curvature for the ensemble of asperities can be obtained from [START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF] by integration

κ g = ∞ 0 ∞ 0 κ 1 κ 2 P(κ 1 , κ 2 ) dκ 1 dκ 2 ≈ 1.356 (12)

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 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 -∆z i + R 2 -(X 0 -x i ) 2 -(Y 0 -y i ) 2 -R -Z 0 = g i0 -∆z i , (13)
where ∆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 + 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 g i 3/2 , a i = r i g i , with c i = 4 3 E * √ r i (14)
where g i = 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:

∆z f = f πE * (x -x f ) 2 + (y -y f ) 2 . ( 15
)
This coupling between additional displacement and forces makes the system of equations ( 13)-( 15) strongly coupled. Then, following the superposition principle, at a given asperity i, its total vertical shift is given by:

∆z i = 1 πE * j∈I\i f j d i j = BF,
where d ij = (x ix j ) 2 + (y iy j ) 2 and the symmetric interaction matrix B (containing the inverse distances) and the vector of forces F are given by:

B = 1 πE *                          0 d -1 12 d -1 13 . . . d -1 1N d -1 12 0 d -1 23 . . . d -1 2N d -1 13 d -1 23 0 . . . d -1 2N . . . . . . . . . . . . . . . d -1 1N d -1 N2 d -1 N3 . . . 0                          , F =                        f 1 f 2 . . . f N                       
Therefore from ( 14) we obtain a non-linear system of equations:

F -C G 0 -BF 3/2 = 0, (16)
where

C = diag[c 1 c 2 . . . c N ]
is the diagonal matrix, and the vector G 0 = [g 10 g 20 . . . g N0 ] 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:

V = [ v 1 v 2 . . . v N ] , V n = [v n 1 v n 2 . . . v n N ]
. 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 x and non-coercive properties. Therefore we regularize it as follows

m(x, δ) =        δ 2 2δ -x , if x < δ x, otherwise,
If δ → 0 we recover the classical Macaulay brackets: m(x, δ) → x . This approximation creates small repulsive forces between the sphere and asperities, but they decay quite rapidly (as z -3/2 ), the parameter δ reduces to zero during the convergence loop enabling to ensure a better convergence for Newton-Raphson method, initial value of δ is chosen to be σ/100. With this smoothing the system of equations ( 16) reduces to

F -Cm(G 0 -BF, δ) 3/2 = 0, (17)
whose solution tends to the solution of ( 16) when δ → 0. The initial guess for the force vector F 0 is computed from a model without interaction (B = 0), then from Eq. ( 16):

F 0 = C G 0 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 ∆F as:

∆F = -K -1 R, (18)
where the free term R and the tangent matrix K are given by:

R = F 0 -Cm(G 0 -BF 0 , δ) 3/2 , K = I - 3 2 Cm(G 0 -BF 0 , δ) 1 2 ∂m(G 0 -BF, δ) ∂F F=F 0
, where I is identity matrix. After every iteration, according to the classical Newton's method, the vector of forces is updated as F i = F i-1 + ∆F and all negative forces are removed as F i = F i . Regardless the fact that function m(x, δ) 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

∆F i / F i ≤ ε,
the tolerance in all simulations is chosen to be ε = 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 ν = 0.42 and the Young's modulus is E = 78 GPa, giving the effective modulus E * ≈ 94.71 GPa. The size of the square region on elastic half-space with identified asperities is set to be 1 mm × 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 separation5 u z = 4σ down to u z = σ 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. 2 , and (2) which are located in the circle (x -X 0 ) 2 + (y -Y 0 ) 2 < r 0 , where r 0 = 1.5 2(4σ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σu z ).

z i > u z + R -R 2 -(X 0 -x i ) 2 -(Y 0 -y i )
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 × N created over the indentation area L × L with cells of size ∆x × ∆x. In each cell {i, j}, i, j ∈ [0, N -1] all contacting asperities I i j ⊂ I contribute to the coarse-grained contact area fraction A ij , and to the coarse-grained contact pressure p i j , where

A i j = 1 ∆x 2 k∈I ij πa 2 k , p i j = 1 ∆x 2 k∈I ij F k , (19)
where a k is computed using ( 14) and I i j is a set of asperities with in-plane coordinates located in the given cell.

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 [START_REF] Bowden | The Friction and Lubrication of Solids[END_REF], 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 [START_REF] Greenwood | The elastic contact of rough spheres[END_REF], 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 [START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF]. However, as remarked in [START_REF] Greenwood | The contact of two nominally flat rough surfaces[END_REF][START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF], 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. 4). The data were averaged over 120 simulations per each combination of root mean squared roughness σ = {0.001, 0.01, 0.1} and different upper cut-offs L/λ s = {256, 512}, the lower cut-off was kept constant L/λ l = 24. The vertical displacement of the indenter's tip was changed from Z 0 = 4σ to σ 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 ∼ N 2/3 . The second regime is characterized by an approximately linear evolution of the contact area A II ∼ N and it corresponds to multi-asperity statistical 20),( 21)) are also plotted.

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 asperities7 r = rg ≈ 0.737/ √ m 4 (the last approximate equality follows from ( 12)):

A I = π 3r g 4E * 2 3 N 2 3 ≈ π 0.55275 √ m 4 E * 2 3 N 2 3 . ( 20
)
In the second regime, at moderate loads, the contact area should grow accordingly to classical multiasperity theories [START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF][START_REF] Bush | The elastic contact of a rough surface[END_REF][START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load?[END_REF]], Persson's theory [START_REF] Persson | Theory of rubber friction and contact mechanics[END_REF] and more recent full-scale numerical simulations [START_REF] Hyun | Finite-element analysis of contact between elastic self-affine surfaces[END_REF][START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF][START_REF] Putignano | The influence of the statistical properties of self-affine surfaces in elastic contacts: A numerical investigation[END_REF][START_REF] Afferrante | Interacting and coalescing hertzian asperities: A new multiasperity contact model[END_REF] with the scaling depending on the root mean squared roughness gradient |∇z| 2 , which for Gaussian isotropic surface can be estimated through the second spectral moments |∇z| 2 = √ 2m 2 (see A). Therefore, in the second regime the contact area growth is determined by the following equation:

A II = βκ √ 2m 2 E * N, (21)
where the constant β ≈ 0.2 was identified from our simulations to fit approximately the data. The factor κ = √ 2π appears from statistical multi-asperity models [START_REF] Bush | The elastic contact of a rough surface[END_REF][START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF][START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load?[END_REF] and was shown to be very accurate for predicting the contact area evolution at small loads for nominally flat surfaces [START_REF] Yastrebov | The role of the roughness spectral breadth in elastic contact of rough surfaces[END_REF]. For higher loads, both multi-asperity models [START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load?[END_REF][START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF] and numerical simulations [START_REF] Paggi | The coefficient of proportionality κ between real contact area and load, with new asperity models[END_REF][START_REF] Yastrebov | The role of the roughness spectral breadth in elastic contact of rough surfaces[END_REF] demonstrate a slightly non-linear area growth which depends on Nayak's parameter α = m 0 m 4 /m 2 2 : greater α 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 [START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF].

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 π 3 (2m 2 ) 3/2 E * β 3 κ 3 m 4 ≈ 212.66 m 3/2 2 E * m 4 , (22)
where the numerical factor 212.66 ≈ 0.55275 2 (π √ 2/κβ) 3 with κ = √ 2π and β = 0.2. Note that in [START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF], a different expression was obtained8 : N * PR = (9π 3 /4m 4 )( √ 2m 2 /κ) 3 E * , which after evaluating all numerical quantities and assuming κ = √ 2π gives N * PR ≈ 12.5287m 3/2 2 E * /m 4 which is considerably smaller than our value [START_REF] Dieterich | Imaging surface contacts: power law contact distributions and contact stresses in quartz, calcite, glass and acrylic plastic[END_REF]. This is because of slight overestimation made by the authors [START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF] of the asperity radius rg = 2/ √ m 4 and of non-taking into account the scaling factor 9 β 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 λ l affects slightly the numerical values, but keeps the transition interval unchanged. The normalized data plotted in Fig. 5 include the results found for different indenter's radius R = {2.5; 5} mm, and different lower cut-off wavelength λ 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 ∼ m 2 ∼ m 0 = σ 2 then the contact areas in the first [START_REF] Derjaguin | Effect of contact deformations on the adhesion of particles[END_REF] and the second ( 21) regimes scale as follows:

A I (N, σ) = A I (N/σ) ∼ N σ 2/3 , A II (N, σ) = A II (N/σ) ∼ N σ . (23)
The scaled results and theoretical predictions are plotted in Fig. 5. 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 * ∼ σ.

Maps of contact area/pressure

The linear regime (21) 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) = κp(x, y)dA 0 /( √ 2m 2 E * ) and the pressure is distributed according to Hertz theory p(x, y) = p 0 1 -(x 2 + y 2 )/a 2 , with p 0 = 3N/(2πa 2 0 ), 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 th = κ √ 2m 2 E * N. (24)
However, the contact area for indentation of the rough contact [START_REF] Greenwood | The elastic contact of rough spheres[END_REF] can be significantly larger than the one predicted by Hertzian theory in case of light indentations, therefore the separation between the Normalized force, N / √2m 2 (N) True contact area, A (mm 2 ) σ = 0.001 mm, σ = 0.01 mm, σ = 0.01 mm, σ = 0.1 mm,

N 2/3 N 2/3 N I ~N 2/3
A II ~N σ = 0.01 mm, R = 5 mm, = 0.01 mm, L/λ l = 48

Figure 5: 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 σ = {0.001, 0.01, 0.1} and different upper cut-offs L/λ s = {256, 512}; if it is not stated differently L/λ l = 24 and R = 2.5 mm; to guide the eye, the 2/3-power law and a linear evolution with relevant factors (Eqs. ( 20),( 21)) are also plotted.

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 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 (24) and the numerical results [START_REF] Dick | Fretting modelling with a crystal plasticity model of ti6al4v[END_REF], 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. 6 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, σ = 10 µm, L/λ l = 512, L/λ s = 24, H = 0.8 at different indentations depth ranging from Z 0 /σ = 4 to Z 0 /σ ≈ 1.4, the latter corresponds to the contact force N ≈ 116.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.

To compute the radial distribution of the contact area fraction and contact pressure at different loads, we use the coarse grained data [START_REF] Dapp | Fluid leakage near the percolation threshold[END_REF]: A i j = A (x i , y j ), p i j = p(x i , y j ), where i, j ∈ [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 µm. To find the average radial data, we first interpolated A i j and p i j in between grid points using bi-cubic Bézier interpolation to obtain à (x, y), 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: The resulting radial pressure distributions at different load levels are presented in Fig. 7(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 [START_REF] Greenwood | The elastic contact of rough spheres[END_REF] 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 ≈ 116.6 (N). Therefore the use of the asperity based model can be justified.

A (r) = 1 2π 2π 0 Ã (X 0 + r cos(φ), Y 0 + r sin(φ)) dφ,

Probability distributions

The deformation regime of individual asperities can be determined by the mean contact pressure at every asperity pk = F k /(πa 2 k ). Therefore it could be useful to construct the probability density function of mean contact pressures P( p), which are experienced at individual asperities. To do so, we carried out simulations for different indentation depths Z 0 /σ = {3, 2.5, 2} with R = 2.5 mm, σ = 10 µm, L/λ l = 512, L/λ s = 24, H = 0.8; for every distribution, the data of 2000 simulations were used, which results in ≈ 320 000, 157 000, 39 000 data points for Z 0 /σ = {2, 2.5, 3} indentation depths, respectively. The results for Z 0 /σ = 2 are shown in Fig. 8, the data look the same for other indentation depths. The numerical data fit very accurately the Rayleigh distribution:

P R ( p) = p ∆ 2 exp - p2 2∆ 2 , ( 25
)
where ∆ 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 ≈ 3σ Y , where σ Y is the plastic yield stress [START_REF] Tabor | The hardness of metals[END_REF][START_REF] Johnson | Contact mechanics[END_REF][START_REF] Hill | A theoretical study of the brinell hardness test[END_REF][START_REF] Mesarovic | Spherical indentation of elastic-plastic solids[END_REF][START_REF] Kogut | Elastic-plastic contact analysis of a sphere and a rigid flat[END_REF]. 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 [START_REF] Swadener | The correlation of the indentation size effect measured with indenters of various shapes[END_REF]:

H = H 0 1 + R * R ,
where R * = r/(ρ s b), and r is the Nye factor10 , ρ 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/ √ 2 ≈ 0.288 nm, the density of statistically stored dislocations in strongly deformed region under the indentation could be estimated to be ρ s ≈ 10 16 m -2 , taking r = 2, we obtain R * ≈ 0.69 µm, which is in good agreement with the data presented in [START_REF] Kim | Indentation size effect for spherical nanoindentation on nanoporous gold[END_REF]. Assuming that the mean asperity radius is rg ≈ 0.737/ √ m4 and assuming macroscopic hardness H 0 = 250 MPa, we obtain the micro-asperity hardness to be H = 250 1 + 0.69/r g with rg measured in micrometers. For the given example with L/λ l = 24, L/λ s = 512 and H = 0.8 and L = 1 mm, we obtain rg ≈ 0.18 µm, raising the micro-hardness to H ≈ 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. 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) = 2p/p 2 0 , if p ≤ p 0 0, otherwise, (26)
where p 0 is the maximal pressure at this asperity, which in turn is given by p 0 = 1.5 p. The resulting full pressure distribution is shown in Fig. 9. 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 [START_REF] Goryacheva | Mechanics of discrete contact[END_REF]. 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) = γp α exp -(p/δ) β , (27)
where we identified α ≈ 2/3 and β ≈ 3/2. The fact that α ≈ 2/3 is consistent with recent numerical findings [START_REF] Wang | Gauging persson theory on adhesion[END_REF] where the author identified α ≈ 0.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. The distribution of individual contact spot sizes associated with contacting asperities is of big importance for adhesive wear [START_REF] Bowden | The Friction and Lubrication of Solids[END_REF][START_REF] Frérot | A mechanistic understanding of the wear coefficient: From single to multiple asperities contact[END_REF] as well as for electric and thermal conductivity [START_REF] Greenwood | Constriction resistance and the real area of contact[END_REF][START_REF] Slade | Electrical contacts: principles and applications[END_REF]. The distribution of areas evaluated for same parameters and for indentation depth Z 0 /σ = 2 is plotted in Fig. 10 and follows exponential decay P(A) ∼ exp(-A/a). Our results share the same features with experimental observations [START_REF] Dieterich | Imaging surface contacts: power law contact distributions and contact stresses in quartz, calcite, glass and acrylic plastic[END_REF] and numerical results [START_REF] Hyun | Elastic contact between rough surfaces: Effect of roughness at large and small wavelengths[END_REF][START_REF] Dapp | Meeting the contact-mechanics challenge[END_REF][START_REF] Frérot | A mechanistic understanding of the wear coefficient: From single to multiple asperities contact[END_REF], 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 ≈ 0.5 % of contact area fraction obtained at these loading conditions with our model.

Discussion

The deterministic multi-asperity model inspired from [START_REF] Afferrante | Interacting and coalescing hertzian asperities: a new multiasperity contact model[END_REF] enables to solve precisely the mechanical normal contact of two spheres in presence of roughness. The semi-analytical solution [START_REF] Greenwood | The elastic contact of rough spheres[END_REF] 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 [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF]. 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 [START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF]. 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 [START_REF] Greenwood | The contact of two nominally flat rough surfaces[END_REF], 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 [START_REF] Greenwood | The elastic contact of rough spheres[END_REF][START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF]. 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 [START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF], 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 [START_REF] Derjaguin | Effect of contact deformations on the adhesion of particles[END_REF], JKR [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF] or Maugis-Dugdale model [START_REF] Maugis | Adhesion of spheres: the jkr-dmt transition using a dugdale model[END_REF][START_REF] Johnson | An adhesion map for the contact of elastic spheres[END_REF]. 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 [START_REF] Greenwood | A note on Nayak's third paper[END_REF][START_REF] Afferrante | Interacting and coalescing hertzian asperities: a new multiasperity contact model[END_REF] 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 -(p/p 0 ) 3/2 , 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 [START_REF] Wang | Gauging persson theory on adhesion[END_REF].

The density of asperities can be found as [START_REF] Nayak | Random process model of rough surfaces[END_REF]

D = √ 3 18π m 4 m 2 ,
which by substituting expression for spectral moments [START_REF] Greenwood | Formulas for moderately elliptical hertzian contacts[END_REF] takes the following form [START_REF] Yastrebov | From infinitesimal to full contact between rough surfaces: evolution of the contact area[END_REF]:

D = √ 3k 2 l 24π 1 -ξ 6 3 + ζ 4-2H -1 2 -H 1 -ξ 4 2 + ζ 2-2H -1 1 -H , (31)
which for a surface without a plateau in the spectral density takes a simpler form

D ≈ √ 3k 2 l 24π 1 -H 2 -H ζ 4-2H ζ 2-2H -1 . ( 32
)
Note that in the numerator we readily approximated (ζ 4-2H -1) by ζ 4-2H ; however, this approximation cannot be made in the denominator if H is reasonably high. In Fig. 11 we compare this analytical expression [START_REF] Greenwood | A note on Nayak's third paper[END_REF] with numerically computed asperity densities for surfaces with different Hurst exponents H ∈ [0.1, 0.9] and different magnifications ζ = k s /k l ∈ [3.125, 125], to keep the surfaces representative we fixed the longest wavelength to L/λ l = 8 and considered surfaces without a plateau. Note that to ensure accuracy of asperity detection, for L/λ ≤ 500 (ζ ≤ 62.5) we used surfaces with discretization N = 4096, and for 500 < L/λ s ≤ 1000 (62.5 < ζ ≤ 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/λ l = 8, L/λ s = 1000, reached ≈ 8% for H = 0.9 and raised to ≈ 15% for H = 0.1. The average distance between closest asperities for a surface without plateau is given by:

d = 1/ √ D ≈ 24π √ 3k 2 l (2 -H) (1 -H) (ζ 2-2H -1) ζ 4-2H = λ l 6 √ 3π (2 -H) (1 -H) (ζ 2-2H -1) ζ 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 = d /2. Moreover, at average contact radius a c = d /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 = λ l √ 3 8π (2 -H) (1 -H) (ζ 2-2H -1) ζ 4-2H . (33)

Figure 1 :

 1 Figure 1: 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.

Figure 2 :

 2 Figure 2: 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/∆x = 2048, lower cut-off L/λ l = 4, and root mean squared height σ/L = 0.025: (a) L/λ s = 16, (b) L/λ s = 32, (c) L/λ 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.

Figure 3 :

 3 Figure 3: Comparison of the joint probability density computed numerically (on the left) and analytically, following Greenwood's equation (11); principal curvatures κ 1 , κ 2 are normalized by the average algebraic curvature κa .

N 2/ 3 N 2 Figure 4 :

 324 Figure 4: Evolution of the contact area for different root mean squared roughness σ = {0.001, 0.01, 0.1} and different upper cut-offs L/λ s = {256, 512}; to guide the eye, the 2/3-power law and a linear evolution with relevant factors (Eqs. (20),(21)) are also plotted.

Figure 6 :

 6 Figure 6: 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 /σ ≈ {4, 3.57, 3.14, 2.71, 2.29, 1.86}, respectively. The black square is of size 1 × 1 mm. The white dashed line denotes the outline of the Hertzian contact area's prediction for the mean corresponding contact force N.

Figure 7 :

 7 Figure 7: 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.

Figure 8 :

 8 Figure 8: Mean asperity pressure probability p evaluated over 2000 simulations for Z 0 /σ = 2 shown in normal and log-scale (see inset). Black lines correspond to fitted Rayleigh distribution (25).

Figure 9 :Figure 10 :

 910 Figure 9: Probability density of the total pressure distribution evaluated over 2000 simulations at Z 0 /σ = 2. Black lines correspond to Eq. (27) with α = 2/3 and β = 3/2.

  n i f i c a t i o n , ζ H u r s t e x p o n e n t , H Eq. (A.5) Numerical

Figure 11 :

 11 Figure 11: Comparison of Eq. (32) with numerically evaluated asperity density of generated rough surfaces with L/λ l = 8. Cross points of the analytical solution correspond to the same parameters used in numerical estimation.

We recall that in the vicinity of the tip, a sphere and a paraboloid of revolution of the same curvature are indistinguishable.

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

l.u. is length units.

In the original paper[START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF], the probability density should be divided by a factor of two to fulfil the normalization.

By separation we understand the distance between the sphere's tip and the mean line of the rough surface.

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.

A more accurate estimation only for highest asperities could be obtained using Eq. (11) from[START_REF] Greenwood | A simplified elliptic model of rough surface contact[END_REF] and by correcting it with a factor 1/2.

Lower index "PR" was used to indicate the authors (Pastewka & Robbins) of the mentioned paper[START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF], we show here the combination of Eq.(5) and (6) from this paper.

This scaling factor is implicitly contained in authors'[START_REF] Pastewka | Contact area of rough spheres: large scale simulations and simple scaling laws[END_REF] estimation of the mean contact pressure N/A 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.

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

Acknowledgment

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

Appendix

A Spectral density and asperity density

Consider a surface spectrum Φ with a plateau which corresponds to an isotropic and periodic rough surface with a period L. Following the construction technique [START_REF] Hu | Simulation of 3-D random rough surface by 2-D digital filter and fourier analysis[END_REF] for all possible directions determined by wavenumbers k x , k y in OX and OY directions, respectively, in average for K = k 2

x + k 2 y , the following equation holds:

where

where from it follows ξ ≤ 1, ζ ≥ 1. A spectral moment m p is given by the following convolution

where T(p) T(p) =