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 