Roughness-Induced Adhesive Hysteresis in Self-Afﬁne Fractal Surfaces

: It is known that in the presence of surface roughness, adhesion can lead to distinct paths of loading and unloading for the area–load and penetration–load relationships, thus causing hysteretic loss. Here, we investigate the effects that the surface roughness parameters have on such adhesive hysteresis loss. We focus on the frictionless normal contact between soft elastic bodies and, for this reason, we model adhesion according to Johnson, Kendall, and Roberts (JKR) theory. Hysteretic energy loss is found to increase linearly with the true area of contact, while the detachment force is negligibly inﬂuenced by the maximum applied load reached at the end of the loading phase. Moreover, for the micrometric roughness amplitude h rms considered in the present work, adhesion hysteresis is found to be affected by the shorter wavelengths of roughness. Speciﬁcally, hysteresis losses decrease with increasing fractal dimension and cut-off frequency of the roughness spectrum. However, we stress that a different behavior could occur in other ranges of roughness amplitude.


Introduction
The hysteretic dissipation is given by the difference between the work needed to bring two bodies into contact and that required to detach them. Its origin may be related to various phenomena occurring at the contact interface. The main causes of hysteresis are viscoelasticity [1,2], plasticity [3], adhesive elastic instabilities at jump-in and jump-out of contact [4], and surface roughness [5]. In particular, all natural and artificial surfaces are rough at some scale. Therefore, hysteretic losses may affect several technological applications. For example, biomedical devices [6] and structural adhesives [7] must safely adhere to surfaces during their application, but they should be easy to remove for reuse. Moreover, a recent challenge in soft robotics is to create climbing robots with reversible adhesion skills [8].
In contact experiments on soft matter, velocity-dependent dissipations are usually measured during detachment as a consequence of bulk viscoelasticity [9]. In a recent work [5], Dalvi et al. carried out loading-unloading contact experiments between smooth silicone hemispheres and rough nanodiamond substrates. Their experiments were performed at very low velocities (60 nm/s) both for the approach and detachment. Such choice allows the avoidance of velocity-dependent dissipations. However, great adhesion hysteresis was still observed due to the roughness-induced increase in the true contact area. This effect is expected to occur in compliant materials with small root mean square (rms) roughness amplitude (h rms 1 nm) [1], when they are bring in complete contact. Moving from the assumption of full-contact conditions, Dalvi et al. applied Persson and Tosatti (PT) adhesion theory [10] for predicting the magnitude of adhesion hysteresis. They found that the hysteretic dissipation increases almost linearly with the true contact area A, and it is equal to the product between A and the intrinsic surface energy ∆γ, which depends on the interfacial adhesive properties of contacting bodies.
In [11], it is experimentally shown that h rms can both increase and decrease the effective adhesive surface energy with respect to the smooth case. Moreover, numerical simulations of continuum adhesive contacts [12] have shown that there is an optimal h rms that leads to a maximization of the hysteretic loss and pull-off force. Such value of h rms is found when the contact region turns from being simply connected to being multiply connected. Similarly, in [13] it is shown that the effective surface energy ∆γ eff reaches a maximum for a certain h rms that, for vanishing applied pressure, is quite close to the value above which the effective contact area A becomes smaller than the nominal one A 0 . Moreover, the enhancement in the adhesion for small h rms is much larger for H < 0.5, where H is the Hurst exponent, as the roughness-induced increase in the surface area is smaller when For RMS roughness amplitudes of the order of few microns, the true area of contact is expected to be predominantly multiply connected. In such case, partial contact conditions occur and surface roughness leads to a reduction in the true area of contact. This in turn destroys adhesion. However, Kesari et al. [11] found that adhesion hysteresis can also be measured for relative large h rms . Inspired by the experimental findings in [11], Deng and Kesari (DK) [14] developed an analytical model for estimating hysteresis losses in the adhesive elastic contacts under the assumption of large roughness. DK's model captures the increase of adhesion hysteresis with the penetration, which is usually called depthdependent hysteresis. Moreover, in this case, a linear increase of the hysteretic dissipation with the area of contact is observed.
Carbone et al. [15] developed a numerical code based on a Boundary Element Method (BEM) for predicting loading-unloading hysteresis loops in the adhesive elastic contact of fractal self-affine 1D rough profiles. Their simulations were conducted under partial contact conditions, with A/A 0 ranging from 0.25 up to 0.5. Due to adhesion hysteresis, two distinct paths were obtained for loading and unloading curves of the area vs. load relation. In particular, they found two sources of energy dissipation, one occurring at small scales and the second one at large scales.
In [16,17], similar multiasperity models have been developed to estimate adhesion hysteresis. They moved from the pioneering Greenwood and Williamson (GW) model, in which roughness is described by a distribution of identical spherical asperities. Adhesion is then implemented according to the classical theory of Johnson, Kendall, and Roberts (JKR) [18]. In a loading-unloading cycle, each asperity exhibits a hysteretic dissipation, which is due to jump-in and jump-off contact instabilities. The total adhesion hysteresis is returned by the contribution of each asperity. However, such models are based on a simplistic description of the surface roughness and do not take into account the elastic coupling between contact regions.
In this work, we propose an investigation of the adhesive elastic contact of rough surfaces, described by self-affine fractal geometries, with an advanced multiasperity model taking into account lateral interactions of asperities according to the authors of [19,20] and adhesion according to JKR theory. Moreover, the model takes also account of the jump-in and jump-off contact instabilities occurring on each asperity.

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 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), where µ = ∆γ 2 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 (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).

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 i of the elastic half-space at the location of the asperity i is whereŵ i is the displacement due to the elastic interaction between the asperities in contact and is given by [24] where n ac is the number of contact spots and r ij 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 i is done by inverting the JKR relation (2). Then, after a further increment of the approach ∆δ i = z i − w i , being z i the height of the asperity i, the contact radius is increased by the quantity 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ū is computed asū 0 − δ, whereū 0 and δ are the initial separation and the total approach, respectively.

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 and zero otherwise. We have denoted with q L = 2π/L and q 1 = 2π/λ 1 the short and long frequencies cut-off, respectively. The quantity L represents the lateral size of the domain (in this case L x = L y ). Finally, H is the so-called Hurst exponent, which is related to the fractal dimension D f = 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 1 = ζ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. For adhesiveless rough contacts, the area-load relation is known to be linear [27][28][29]. However, recent studies confirm that adhesion may lead to strong non-linearity of the F − A curve [30][31][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].

Adhesive Hysteresis and Pull-Off Force: Effect of Loading Parameters
In an approach-retraction cycle, the magnitude Θ of hysteretic losses is equal to the area between the loading and unloading F − ∆ curves, where ∆ is the mean penetration of the rough surface in the half-space. Figure 3C shows the evolution of the dimensionless contact force F/(E * 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 * A 0 h rms ) with A/A 0 corresponding to the maximum applied loads. The model captures the linear relation between Θ and A/A 0 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.

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 1 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 f = 3 − H). Figure 4A shows the F/(E * A 0 ) − A/A 0 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 * 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 * 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 0 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]. 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 * 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 * A 0 ) − ∆/h rms relation modifies with H. In particular, at H = 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.

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 ζ 10 7 , 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.  For a rough surface the energy loss of each asperity depends on the value of its radius of curvature. Following the work in [43], the variation of curvature in the population of all asperities contained in any unit region is where t = − √ 3/m 4 k m , being k m ∈ (0, ∞) the surface's mean curvature at the apex of an asperity, and the constants C 1 = α/(2α − 3) and C 2 = C 1 √ 12/α are related to the Nayak parameter α = m 0 m 4 /m 2 2 . The quantities m 0 , m 2 , and m 4 are the spectral moments of surface roughness PSD. In particular, they are related to the rms roughness amplitude h rms , gradient h rms and curvature h rms by h rms = √ m 0 , h rms = √ 2m 2 , and h rms = √ 8/3m 4 . Substituting (A2) in (A1) and integrating on the range of variation of t, the mean energy loss of contacting asperities can be written as being R(t) = (1/R tip − √ m 4 /3t) −1 .In DK's model, R tip is the radius of the spherical indenter that is in contact with a nominally flat surface. In our case, as contact occurs between nominally flat surfaces, 1/R tip → 0.
Finally, the total energy loss can be computed as where η = m 4 /(6π √ 3m 2 ) is the asperity density in a nominal contact region of unit area. In the original DK's model, the contact area A is computed applying JKR theory to the macroscopic spherical indenter of radius R tip and neglecting roughness contribution. In particular, A is calculated as A ∆ max − A ∆in , where A ∆ max and A ∆in are the values of the contact area at the maximum indentation and the macroscopic jump-in instability of the spherical tip, respectively.
In this work, as we are dealing with the contact between two nominally flat surfaces, A ∆in is interpreted as the true contact area of the first few asperities jumping into contact, while A ∆ max is the true area of contact obtained at the end of the loading phase. Moreover, as in our calculations A ∆in A ∆ max , we have assumed A ∆ max − A ∆in ≈ A ∆ max .