**2. Problem Statement**

#### *2.1. Adhesion Hysteresis of Smooth Elastic Spheres*

JKR theory [18] is commonly used to predict adhesion of soft materials. In the case of spherical contact, the fundamental equations of JKR theory give the applied load *F* and penetration *δ* as a function of the contact radius *a*

$$F\_{\perp} = \frac{4}{3} \frac{E^\* a^3}{R} - \sqrt{8\pi E^\* \Delta \gamma a^3} \tag{1}$$

$$
\delta\_{\perp} = \frac{a^2}{R} - \sqrt{\frac{2\pi a \Delta \gamma}{E^\*}} \,, \tag{2}
$$

being *R* the radius of curvature, *E* ∗ the composite elastic modulus of the contacting bodies and ∆*γ* the interface adhesion energy.

Under controlled displacement conditions, JKR model predicts a jump into contact at *δ* = 0. However, Wu [21], investigating the jump-in instability occurring in atomic force microscopy measurements, found that jump-in instability is reached at a critical gap *δ*in. Such effect is due to van der Waals interactions acting between approaching bodies. Wu proposed an empirical formula for the jump-in distance (valid for *µ* ≥ 2),

$$
\delta\_{\rm in} = \left(1 - 2.641 \mu^{3/7}\right) \varepsilon \tag{3}
$$

where *µ* = ∆*γ* <sup>2</sup>*R*/*E* ∗ 2 1/3 /*ǫ* is the so-called Tabor parameter [22] and *ǫ* is the range of attractive forces.

The above equation can be used to modify JKR theory to consider the jump-in critical distance (see in [23]).

Moreover, under displacement controlled conditions and during retraction, JKR theory predicts a jump-off instability at a critical penetration

$$\delta\_{\rm off} = -\left(\frac{27\pi^2 \Delta \gamma^2 R}{64E^{\*2}}\right)^{1/3}.\tag{4}$$

Figure 1 shows the loading-unloading cycle predicted by JKR theory. The yellow area represents the energy loss due to jumping instabilities. For smooth contact, the energy loss is independent on the maximum penetration (or, equivalently, applied force).

**Figure 1.** The load–penetration curve predicted by Johnson, Kendall, and Roberts (JKR) theory. The loading (unloading) path is denoted by green (red) arrows. Positive penetration values correspond to indentation and compressive force. Jump into contact occurs at a penetration *δ*in. The unloading path overlaps the loading one, but jump out of contact occurs at a critical penetration *δ*off. The yellow area denotes the hysteretic energy loss.

#### *2.2. Adhesion Hysteresis of Rough Elastic Surfaces*

In the case of rough surfaces, multiple unstable jumps occur at the location of each asperity in a loading–unloading cycle. To take account of this phenomenon, Equations (3) and (4) are implemented in a multiasperity model which in turn takes into account the elastic coupling due to asperities lateral interactions.

Let us consider a rigid rough surface approaching an elastic half-space (Figure 2). According to JKR formalism, the normal displacement *w<sup>i</sup>* of the elastic half-space at the location of the asperity *i* is

$$w\_i = \frac{a\_i^2}{R\_i} - \sqrt{\frac{2\pi a\_i \Delta \gamma}{E^\*}} + \mathfrak{w}\_i \tag{5}$$

where *w*ˆ*<sup>i</sup>* is the displacement due to the elastic interaction between the asperities in contact and is given by [24]

$$\begin{array}{rcl} \mathfrak{w}\_{i} &=& \sum\_{j=1, j \neq i}^{n\_{\rm ac}} \frac{a\_{j}^{2}}{\pi R\_{j}} \left( \sqrt{\frac{r\_{ij}^{2}}{a\_{j}^{2}}} - 1 + \left( 2 - \frac{r\_{ij}^{2}}{a\_{j}^{2}} \right) \arcsin\left(\frac{a\_{j}}{r\_{ij}}\right) \right) \end{array} \tag{6}$$

$$-\frac{1}{\pi a\_j E^\*} \sqrt{8\pi a\_j^3 E^\* \Delta \gamma} \arcsin\left(\frac{a\_j}{r\_{ij}}\right) \text{ for } r\_{ij} > a\_j \tag{7}$$

$$\left| \mathfrak{w}\_{\dot{l}} \right| = \left| \frac{a\_{\dot{j}}^2}{R\_{\dot{j}}} - \sqrt{\frac{2\pi a\_{\dot{j}}\Delta\gamma}{E^\*}} \right. \tag{8}$$

where *nac* is the number of contact spots and *rij* is the distance between the asperities *i* and *j*.

When the rough surface approaches the half-space, a new contact is formed when the gap between an asperity and the half-space becomes smaller than *δ*in, which is calculated for each asperity. A first estimate of the asperity contact radius *a<sup>i</sup>* is done by inverting the JKR relation (2). Then, after a further increment of the approach <sup>∆</sup>*δ<sup>i</sup>* = *<sup>z</sup><sup>i</sup>* − *<sup>w</sup><sup>i</sup>* , being *z<sup>i</sup>* the height of the asperity *i*, the contact radius is increased by the quantity

$$
\Delta a\_{\rm i} = \frac{\Delta \delta\_{\rm i}}{2a\_{\rm i}/R\_{\rm i} - \sqrt{\pi \Delta \gamma/(E^\* a\_{\rm i})}} \tag{9}
$$

which is obtained by differentiating Equation (2).

The total contact area and load are then obtained by summing up the contributions of all the asperities in contact. Moreover, as a self-balanced load distribution is considered, the interfacial mean separation *u*¯ is computed as *u*¯<sup>0</sup> − *δ*, where *u*¯<sup>0</sup> and *δ* are the initial separation and the total approach, respectively.

**Figure 2.** Elastic half-space in contact with a rigid rough surface.

## **3. Results**

Computations have been performed on fractal self-affine isotropic surfaces. Roughness is described by its power spectral density (PSD), which has a power law relation with the magnitude *q* = |**q**| of the wavevector **q** . In this work, we consider fractal surfaces with PSD

$$\mathcal{C}(q) = \mathcal{C}\_0 (q/q\_\mathcal{L})^{-2(H+1)} \qquad \qquad \text{for } q\_\mathcal{L} \le q < q\_\mathcal{1} \tag{10}$$

and zero otherwise. We have denoted with *q*<sup>L</sup> = 2*π*/*L* and *q*<sup>1</sup> = 2*π*/*λ*<sup>1</sup> the short and long frequencies cut-off, respectively. The quantity *L* represents the lateral size of the domain (in this case *L<sup>x</sup>* = *Ly*). Finally, *H* is the so-called Hurst exponent, which is related to the fractal dimension *D*<sup>f</sup> = 3 − *H*. Rough surfaces are numerically generated according to the spectral method proposed in [25,26].

In our calculations, we fixed *L* = 1 mm and *h*rms = 5 µm. Furthermore, two sets of simulations have been performed. In the first one, we fixed *H* = 0.8 and *q*<sup>1</sup> = *ζq*L, with magnification *ζ* = 64, 128, 256, 512. In the second one, we fixed *ζ* = 128 and *H* = 0.45, 0.65, 0.8, 0.95.

#### *3.1. Adhesive Hysteresis and Pull-Off Force: Effect of Loading Parameters*

Figure 3A shows the normalized true contact area *A*/*A*<sup>0</sup> as a function of the dimensionless load *F*/(*A*0*E* ∗ ). Calculations have been performed for *H* = 0.8, *ζ* = 128, and *h*rms = 5 µm. Moreover, the curves are obtained by averaging the results of six surface realizations. The material properties are *E* <sup>∗</sup> = 1.0 MPa and ∆*γ* = 0.07 J/m<sup>2</sup> , which are typical values for very soft silicon elastomers. Unloading starts from different maximum applied loads *F*max/(*E* <sup>∗</sup>*A*0) = 0.004, 0.0071, 0.012, 0.015.

For adhesiveless rough contacts, the area–load relation is known to be linear [27–29]. However, recent studies confirm that adhesion may lead to strong non-linearity of the *F* − *A* curve [30–32]. In our calculations, this is especially true for the unloading path in agreement with numerical [15] and experimental [33] findings. The pull-off force is the maximum negative load reached during retraction. Near the pull-off point, unloading paths almost collapse on a single curve. As a result, the pull-off force *F*po is quite independent on the maximum true area of contact reached during the approach (Figure 3B). This is in agreement with experimental findings of Refs. [11,33].

In an approach–retraction cycle, the magnitude Θ of hysteretic losses is equal to the area between the loading and unloading *<sup>F</sup>* − <sup>∆</sup> curves, where <sup>∆</sup> is the mean penetration of the rough surface in the half-space. Figure 3C shows the evolution of the dimensionless contact force *F*/(*E* <sup>∗</sup>*A*0) with the normalized penetration ∆/*h*rms. Recent experiments [5,11] suggest that, in presence of surface roughness, Θ increases linearly with the penetration ∆ (or in a similar way with the true area of contact *A*). Such phenomenon is known in literature as depth-dependent hysteresis. Figure 3D shows the increase of the dimensionless energy loss Θ/(*E* <sup>∗</sup>*A*0*h*rms) with *A*/*A*<sup>0</sup> corresponding to the maximum applied loads. The model captures the linear relation between Θ and *A*/*A*<sup>0</sup> and the analytical predictions given by DK's model [14], as adapted to the present case (see Appendix A), are coherent with our numerical results.

**Figure 3.** (**A**) The normalized area of contact *A*/*A*<sup>0</sup> as a function of the dimensionless load *F*/(*E* ∗*A*0). Results are obtained on surfaces with *h*rms = 5 µm, *H* = 0.8 and *ζ* = 128. Loading (green dashed line) and unloading (red solid line) curves are shown. Results are averaged on 6 surface realizations. (**B**) The dimensionless pull-off force *F*po/(*E* ∗*A*0) as a function of the maximum applied load at the end of the loading phase. Error bars denote the standard deviation on 6 surface realizations. (**C**) The dimensionless mean penetration ∆/*h*rms as a function of the applied load *F*/(*E* ∗*A*0). (**D**) The dimensionless energy loss Θ/(*E* ∗*A*0*h*rms) as a function of the normalized area of contact *A*/*A*0. Red solid and blue dashed lines refer to the present calculations and DK's predictions, respectively. Error bars denote the standard deviation on 6 surface realizations.

#### *3.2. Adhesive Hysteresis and Pull-Off Force: Effect of the Fractal Parameters*

Surface roughness can be described by its statistical parameters, i.e., RMS roughness amplitude *h*rms, RMS gradient *h* ′ rms, and RMS curvature *h* ′′ rms. The first one is related to low frequencies of the PSD spectrum, while RMS slope and curvature mainly depends on the cut-off frequency *q*<sup>1</sup> and therefore on the magnification *ζ*. Increasing *ζ*, the PSD spectrum is enriched by smaller and smaller roughness wavelengths. An other important parameter is the Hurst exponent *H*. Low (high) values of *H* correspond to high (low) fractal dimension *D*f (*D*<sup>f</sup> = 3 − *H*).

Figure 4A shows the *F*/(*E* <sup>∗</sup>*A*0) − *A*/*A*<sup>0</sup> relation at increasing values of the magnification *ζ* (*ζ* = 64, 128, 256, 512) for *h*rms = 5 µm and *H* = 0.8. All curves are obtained for a same value of the applied load *F*/(*E* <sup>∗</sup>*A*0) = 0.071 reached at the end of the approach, in similar way to the experiments performed in [5]. The true area of contact decreases with *ζ* as an increase in *ζ* corresponds to bigger rms gradient *h* ′ rms. In such case, surface roughness is described by several length-scales and a greater load is required to create new contact patches on smaller wavelengths. Specifically, as shown in [34], the dependence of the curves on *ζ* (and thus on *h* ′ rms) is exclusively due to the contribution of the repulsive interactions.

Figure 4B shows the dimensionless pull-off force *F*po/(*E* <sup>∗</sup>*A*0) as a function of the magnification. A general drop in the pull-off force is observed by increasing *ζ*. This is in agreement with recent numerical findings in [35], where an in-house Boundary Element Method (BEM) has been developed for studying adhesive contact of rough surfaces, by including full Lennard–Jones potentials and surface integration at the asperity level. However, such results are due to the fact that we are considering surfaces with roughness amplitude of the order of microns. In fact, as found in [13], at sufficiently high values

of the rms roughness amplitude (and more precisely of the product *q*0*h*rms, being *q*<sup>0</sup> the roll-off frequency), for *H* = 0.8 a decrease in magnification *ζ* involves an increase in the effective surface energy at short length scale (large *ζ*). This effect results from the increase in the contact area as more a more short-wavelength roughness components are taken into account. However, we stress that such a result works as long as *q*0*h*rms is high enough. Indeed, at lower *h*rms, such effect is not observed and the adhesion seems to be governed only by the surface roughness amplitude [36].

**Figure 4.** (**A**) The normalized area of contact *A*/*A*<sup>0</sup> as a function of the dimensionless load *F*/(*E* ∗*A*0). Results are obtained on surfaces with *h*rms = 5 µm, *H* = 0.8, and *ζ* = 64, 128, 256, 512. Loading (green dashed line) and unloading (red solid line) curves are shown. Results are averaged on 6 surface realizations. (**B**) The dimensionless pull-off force *F*po/(*E* ∗*A*0) as a function of the magnification. Error bars denote the standard deviation on 6 surface realizations. (**C**) The dimensionless mean penetration ∆/*h*rms as a function of applied load *F*/(*E* <sup>∗</sup>*A*0). (**D**) The dimensionless energy loss Θ/(*E* ∗*A*0*h*rms) as a function of the magnification. Red solid and blue dashed lines refer to the present calculations and DK predictions, respectively. Error bars denote the standard deviation on 6 surface realizations.

In fact, for large *ζ* and small *h*rms, a model based on a JKR-type approach becomes questionable as the dimension of the contact spots decreases [37]. In such case, a DMT-type approach, based on the assumption of long-range adhesion interactions, could be more accurate in modeling the contact problem. In this regard, a DMT-type model is developed in [34] with the aim of investigating the adhesive contact of surfaces with *h*rms of the order of 1 nm. It is found that *F*po is almost independent of *h* ′ rms, i.e., the adhesion force required for the detachment is magnification independent.This is confirmed in [38], where a stickiness criterion [38] is derived from Persson–Scaraggi DMT theory [37]. For typical values of the Hurst exponent (*H* > 0.6), the criterion suggests that adhesion is destroyed by the long wavelengths of roughness, while *ζ* has negligible effects. Such result has been corroborated by very recent experimental [39] and analytical works [40,41], according to which the main parameter "killing" adhesion seems to be the roughness amplitude *h*rms.

Figure 4C shows the load-penetration relationship for increasing magnification *ζ*. The corresponding hysteretic losses Θ are shown in Figure 4D. Our simulations suggest that the magnitude of energy loss follows the same trend of the pull-off force, i.e., it reduces with *ζ*. Once again, DK's predictions, as given by the proposed modified equation in appendix, are in agreement with our numerical calculations. The error bars in Figure 4B–D show the standard deviation on six surface realizations. The scatter is larger for surfaces with low magnification *ζ* as a result of the smaller number of surface asperities. On the

contrary, increasing *ζ*, smaller and smaller asperities are added to the rough surface and their spatial distribution is expected to be more uniform in the nominal contact region. Moreover, as the linear size of the system is finite the PSD is not continuous and even assuming spectral components with random phases uniformly distributed in the range 0 < *φ* < 2*π* the surface will be not ergodic, and a single realization of the surface will be in general highly non-Gaussian, thus entailing finite-size effects related to the finite value of the asperities heights. It follows that for surfaces without a low-wavenumber roll-off (or cut-off) region, quantities which depend on the long wavelength roughness, such as the average interfacial separation (and hence the hysteresis dissipation) at low contact pressures, will vary strongly from one realization to another.

A second set of simulations has been performed on surfaces with fixed *h*rms = 5 µm, *ζ* = 128 and different Hurst exponents *H* = 0.45, 0.65, 0.8, 0.95. We have fixed again the maximum applied load *F*/(*E* <sup>∗</sup>*A*0) = 0.071. An increase in *H* leads to a decrease in RMS gradient, thus explaining why the area increases with *H* for a fixed load (Figure 5A). Figure 5B shows that the pull-off force is destroyed at low *H*; the same trend has been observed in [42], where the adhesive contact between a parabolic indenter with superimposed roughness and an elastic half space has been studied in the JKR-limit. Figure 5C shows how the *F*/(*E* <sup>∗</sup>*A*0) − <sup>∆</sup>/*h*rms relation modifies with *<sup>H</sup>*. In particular, at *<sup>H</sup>* = 0.45 loading and unloading paths overlaps thus showing vanishing adhesive hysteretic loss Θ, which is strongly affected by the Hurst exponent as shown in Figure 5D. Notice results are more scattered for surfaces with lower fractal dimension as finite-size effects related to the absence of a low-wavenumber roll-off (or cut-off) region are exaggerated at higher values of *H*.

**Figure 5.** (**A**) The normalized area of contact *A*/*A*<sup>0</sup> as a function of the dimensionless load *F*/(*E* ∗*A*0). Results are obtained on surfaces with *h*rms = 5 µm, *ζ* = 128, and *H* = 0.45, 0.65, 0.8. Loading (green dashed line) and unloading (red solid line) curves are shown. Results are averaged on 6 surface realizations. (**B**) The dimensionless pull-off force *F*po/(*E* ∗*A*0) as a function of the Hurst exponent. Error bars denote the standard deviation on 6 surface realizations. (**C**) The dimensionless mean penetration ∆/*h*rms as a function of applied load *F*/(*E* <sup>∗</sup>*A*0). (**D**) The dimensionless energy loss Θ/(*E* ∗*A*0*h*rms) as a function of the Hurst exponent. Red solid and blue dashed lines refer to the present calculations and DK predictions, respectively. Error bars denote the standard deviation on 6 surface realizations.

## **4. Discussion and Conclusions**

In the adhesive contact between an elastic half-space and a rigid randomly rough surface, loading-unloading loops can be observed as a result of adhesion hysteresis induced by roughness. Hysteretic losses are found to be linearly increasing with the true area of contact reached at the end of the loading path. On the contrary, the pull-off force is negligibly influenced by the maximum contact area (and thus maximum applied load).

Pull-off force and hysteretic losses are strongly affected by roughness parameters. Specifically, here we have investigated the effects of the Hurst exponent *H* and magnification *ζ*.

Detachment force and hysteretic losses are observed to reduce by decreasing *H* and increasing *ζ*. Such results are related to the increase in the RMS gradient *h* ′ rms occurring when *H* is reduced or *ζ* is increased. Our outcomes are in agreement with the trends shown by very recent numerical, experimental and analytical findings.

However, we stress that our results are obtained on surfaces with RMS roughness amplitude of the order of few micrometers where we can reasonably expect partial contact conditions occur in a wide range of applied loads. In fact, multiasperity models become progressively less accurate moving towards full contact conditions. Moreover, numerical models allow to consider a limited range of magnifications, while real surfaces are characterized by roughness on several length scales (with the modern technologies we can measure *<sup>ζ</sup>* <sup>≃</sup> <sup>10</sup><sup>7</sup> , ranging from centimeter to nanometer scales). Despite such limitations, the present findings help to clarify some aspects of the hysteretic phenomenon occurring in the adhesive contact of rough soft matter.

**Author Contributions:** G.V. and L.A. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the project "FASTire (Foam Airless Spoked Tire): Smart Airless Tyres for Extremely-Low Rolling Resistance and Superior Passengers Comfort" funded by the Italian MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) call 2017—grant n. 2017948FEN.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data presented in this work are available on request from the corresponding author.

**Acknowledgments:** L.A. and G.V. acknowledge support from the Italian Ministry of Education, University and Research (MIUR) under the program "Departments of Excellence" (L.232/2016).

**Conflicts of Interest:** The authors declare no conflicts of interest.
