Finite deformations govern the anisotropic shear-induced area reduction of soft elastic contacts

Solid contacts involving soft materials are important in mechanical engineering or biomechanics. Experimentally, such contacts have been shown to shrink significantly under shear, an effect which is usually explained using adhesion models. Here we show that quantitative agreement with recent experiments can be obtained, with no adjustable parameter, using a non-adhesive model, provided that finite deformations are taken into account. Analysis of the model uncovers the basic mechanisms underlying shear-induced area reduction, local contact lifting being the dominant one. We confirm experimentally the relevance of all those mechanisms, by tracking the shear-induced evolution of tracers inserted close to the surface of a smooth elastomer sphere in contact with a smooth glass plate. Our results suggest that finite deformations are an alternative to adhesion, when interpreting a variety of sheared contact experiments involving soft materials.


Introduction
Rough contacts are ubiquitous in both natural and engineering systems, and indeed, rough contact mechanics have been actively investigated in the last decades (see e.g. Vakis et al. (2018) for a recent review). Most of the effort has been devoted to the normal contact of frictionless interfaces (see Müser et al. (2017) for a comparison of various modelling approaches to such a problem). However, recent experiments involving soft materials like polymers or human skin have revealed complex changes to the contact morphology when a frictional rough contact is submitted to an additional shear load. Not only is the overall real contact area significantly reduced (Sahli et al., 2018;Weber et al., 2019), but it also becomes increasingly anisotropic , two effects that have not been satisfactorily explained yet. For both effects, at the contact interface between linear elastic solids are uncoupled (Johnson, 1985), so that the observed effect of the tangential load on the contact area is unexpected.
Here, we investigate the hypothesis that the shear-induced contact area reduction is an elastic effect enabled by the nonlinear, finite-deformation behaviour of the elastomer. Recently, this possibility was suggested in Wang et al. (2020), and some preliminary support was brought by Mergel et al. (2020) using 2D simulations, but it has never been tested on 3D sphere/plane contacts. To fully test the hypothesis, we have developed a computational 3D model that combines a neo-Hookean hyperelastic model of the sample and a non-adhesive contact model in which friction is governed by a (regularized) Tresca model. This computational model suitably fits the macroscopic experimental data from Sahli et al. (2018Sahli et al. ( , 2019 and quantitatively reproduces the anisotropic evolution of contact area, with no adjustable parameter (Section 2). The presented model results strongly suggest that the contact area reduction is governed by the finitedeformation mechanics, with a dominating effect of local contact lifting and less pronounced effect of micro-slip-induced in-plane deformation. To test those model conclusions, we then present original experiments in which those two mechanisms can indeed be directly observed and quantified (Section 3).

Finite-strain modelling of shear-induced contact area reduction
The finite-strain framework and the Tresca friction model are the two essential features of our model, and these are combined with the finite-element method as a suitable spatial discretization scheme. Since the individual ingredients of the model are rather standard, perhaps except for some details of the computational treatment, the model is only briefly summarizeded in Section 2.1, and its detailed description is provided in Appendix A.

Computational model: finite-strain framework and Tresca friction
The contact problem under consideration, sketched in Fig. 1, corresponds to the normal and tangential loading of a hyperelastic spherical solid in (adhesionless) frictional unilateral contact with a rigid plate. The finite-strain framework employs the geometrically exact kinematics of finite deformations (Appendix A.1) and contact (Appendix A.2), as well as adequate constitutive descriptions of both the bulk elasticity of the sphere (Appendix A.1) and the frictional behaviour at the contact interface (Appendix A.3).
Hyperelasticity is treated using the Mooney-Rivlin model, in the nearly incompressible version. Disregarding compressibility and the related bulk modulus, the model (Eq. (A.3)) involves two parameters, µ 1 and µ 2 , such that µ = µ 1 + µ 2 , with µ the shear modulus. In the following, we will consider two particular cases reducing to a single parameter, µ: (i) the neo-Hookean model, for which µ 2 = 0 and µ 1 = µ, and (ii) an arbitrarily chosen other combination of parameters, namely µ 1 = µ 2 = µ/2 (simply denoted as Mooney-Rivlin from now). Note that µ can be readily determined for any particular experimental sphere/plane contact from the knowledge of the contact area, A 0 = A(Q = 0), under a pure, known normal load, P.
For a wide range of tribological material pairs, the static friction force of sphere/plane contacts is measured to be proportional to the contact area (see Sahli et al. (2018) and references therein). It is thus appealing to examine the simplest friction model that automatically produces such a dependence, namely the Tresca friction model, which is used in our model (Appendix A.3). In the Tresca model, no slip occurs locally until the contact shear strength, σ, is reached, rigid plate rigid support hyperelastic sample (a) trailing edge leading edge (b) Figure 1: 2D sketch of the (3D) soft-sphere/rigid-plane contact under study. The initial configuration is indicated by dashed lines, the deformed configuration is denoted by solid lines for two stages: (a) after the normal load P is applied, and (b) when an additional tangential displacement d is applied, giving rise to a tangential load Q. and σ is assumed constant and independent on the normal contact traction. 1 The Tresca model, even if oversimplified, is expected to provide a reasonable approximation of the tangential stress distribution (at least) when the static friction force, Q s , is reached. Note that σ can be measured in all experiments that monitor the contact area (e.g., Savkoor and Briggs, 1977;Waters and Guduru, 2010;Sahli et al., 2018;Mergel et al., 2019), as the ratio of Q s over A s = A(Q = Q s ).
The finite-element treatment of the contact problem at hand is described in Appendix A.4. Note that the Tresca friction model may lead to significant convergence problems, particularly when the interfacial shear strains are large, which will be the case in our simulations. To circumvent those problems, a Prakash-Clifton-like regularization scheme (Appendix A.3) has been employed, which allowed our actual computations to be successfully carried out for a relatively fine finite-element mesh.
We have applied our model to perform a direct quantitative comparison with the PDMSsphere/glass-plate experiments reported in Sahli et al. (2018Sahli et al. ( , 2019. In particular, the geometry, boundary conditions and loading conditions match the experimental ones. The elastomer sample (identical to that in inset of Fig. 9) is a cylinder of diameter 30 mm and thickness 6 mm, the top of which features a 2 mm-thick spherical cap with a radius of curvature of 9.42 mm. It is fixed on a rigid support. The spherical cap is first brought into normal contact with a rigid plate under constant normal load P (we chose four representative values of P covering the whole range explored in Sahli et al. (2018): 0.27, 0.55, 1.65 and 2.12 N). The contact is then sheared by pulling the plate horizontally with a constant velocity V = 0.1 mm/s. A value of µ = 0.60 MPa was derived in a unique manner using the experimental values of P and A(Q = 0) from the initial contact (see Appendix A.5). Similarly, a value of σ = 0.41 MPa was derived using the experimental values of A s and Q s . Note that the ratio σ/µ is significant, such that shear strains exceeding 50% are expected. This is far beyond the range of validity of the small-strain theory, and linear elasticity in particular, so that the finite-strain framework used here is actually essential. Figure 2 shows the finite-element mesh used in the case of P = 2.12 N (note that the size of the refined-mesh region is adjusted to the normal load), and zoom-ins for select configurations. In particular, panel (d) illustrates the deformation pattern for Q = Q s , showing a significant shear deformation in the subsurface layer and a non-uniform in-plane deformation of the contact surface. The corresponding high mesh distortion, visible at the leading edge, and also shown in detail in Fig. A.17(a), is due to the jump of the tangential contact traction, introduced by the Tresca model, and due to the non-conforming discretization of the contact area boundary (see the related discussion in Appendix A.5).

Shear-induced contact area reduction
Using the computational model described above, we simulated the shear-induced contact area reduction in the sphere/plane experiments of Sahli et al. (2018Sahli et al. ( , 2019, for the two hyperelastic models (neo-Hookean and Mooney-Rivlin). Typical results are reported in Fig. 3(a) which shows the concurrent evolutions of the tangential force, Q, and the contact area, A, as a function of the displacement of the rigid plate, d. For comparison, the results corresponding to linear elasticity are also included in Fig. 3, even if the range of strains expected in our conditions does not really admit application of the small-strain framework. Figure 3(a) shows that the three models yield drastically different predictions, both for the tangential force and the contact area. While the initial contact area at zero tangential force is hardly affected by the model, major differences appear as the tangential force increases. A significant area reduction is predicted for both hyperelastic models, but it remains negligible for the linear elastic model. This result already indicates that the finite-deformation framework is actually a prerequisite to reproduce an area reduction of several tens of percents, as observed experimentally. Note that the amplitude of the reduction is larger for the Mooney-Rivlin model (37% for P = 0.27 N) than for the neo-Hookean model (28%), showing that the effect of the material model that governs the hyperelastic behaviour of the sample is significant. Those very different final contact areas, combined with the very same contact shear strength, σ, explain why the final friction force is significantly larger for the linear elastic model, and is also different for both hyperelastic models. 5 : Shear-induced contact area reduction predicted for various elastic models. (a) Concurrent evolution of the tangential force, Q (dots), and the contact area, A (squares), as a function of the rigid plate displacement, d, for three different elastic models: linear elastic (dotted blue), neo-Hookean (solid purple) and Mooney-Rivlin (dashed green). P = 0.27 N. (b) A vs Q, for the three elastic models and for the four normal loads. The markers correspond to selected instants; the computations were carried out with a much finer time stepping. The solid straight line in panel (b) goes through the origin and has a slope of 1/σ. Figure 3(b) synthesizes all results by showing, for all normal loads and all three elastic models, the contact area A as a function of the tangential force Q. The absence of area reduction for the linear elastic model and the difference in its amplitude for the two hyperelastic models is clearly evidenced. The fact that all curves end on the red line indicates that a homogeneous shear stress distribution equal to the contact shear strength σ is suitably enforced at the onset of gross sliding. Recalling that the material parameters for the Mooney-Rivlin model (µ 1 = µ 2 = µ/2) have been adopted here quite arbitrarily, and that the neo-Hookean is another particular case (µ 1 = µ, µ 2 = 0), one can expect that a range of responses (between and beyond the two tested models) could actually be obtained by varying µ 1 and µ 2 under the constraints µ 1 +µ 2 = µ, µ 1 > 0 and µ 2 ≥ 0.
The differences between the various elastic models concern also the shape of the contact zone. In Fig. 4(b)-(d), typical model predictions for the final contact shape are shown for the three elastic models (red regions). Although all models start with an almost identical circular contact when Q = 0 (dashed lines), the final contact shape has significant differences. The linear elastic model predicts a final contact which remains circular, with almost the same area as in the initial configuration, consistently with the results of Fig. 3. The large area reduction in the Mooney-Rivlin model occurs while keeping an essentially circular contact shape. In contrast, the neo-Hookean model leads to an ellipse-like contact shape, with the size reduction occurring almost exclusively along the shear direction. The latter observation, combined with the fact that 6 0.5 mm the final size along shear is very similar for both hyperelastic models, qualitatively explains why the Mooney-Rivlin model predicts a significantly larger amplitude of area reduction (Fig. 3).
We are now in a position to compare quantitatively our model results to the experimental results of Sahli et al. (2018Sahli et al. ( , 2019. In Fig. 4, panel (a) shows the experimental counterpart of panels (b)-(d). It clearly appears that the model that most closely matches the final shape in the experiments is the neo-Hookean model. In particular, both the size and eccentricity of the ellipselike contact shape is very-well captured. The excellent agreement of the neo-Hookean model with experiments is further demonstrated in Fig. 5. Figure 5 Fig. 3b). In both cases, the amplitudes and shapes of the curves are well captured, although the model slightly underestimates A and does not capture the slight increase of ⊥ observed in the experiments. Note that a similarly good agreement has been obtained for an alternative model with different regularizations and numerical implementations (see Appendix A.6), thus showing the robustness of our results.

Elementary mechanisms of contact area reduction
In Section 2.2, we have shown that our neo-Hookean model provides, without any adjustable parameter, a very good quantitative prediction of the shear-induced contact area reduction observed in experiments. On this model only, we will now perform a thorough analysis of the simulation results to understand what are the elementary mechanisms responsible for such a reduction.
Important qualitative insights can already be obtained from a careful inspection of Fig. 4. First, the green regions around the final contact zone correspond to all nodes that were initially in contact, and which are not anymore at the onset of macroscopic sliding. Those nodes have been lifted out of contact during the incipient loading phase, due to the shear-induced deformations of the elastic sphere. This is the first elementary mechanism potentially causing area variations (here a decrease). Below, it will be denoted as "contact lifting" (or simply "lifting"). Second, small parts (hardly visible in Fig. 4(b)) of the final contact are found outside the initial contact region (shown as a green region in the current configuration). This may be interpreted as nodes Area of contact, A that were initially out of contact, but came into contact during the incipient loading phase. This is the second elementary mechanism potentially causing area variations (here an increase). Below, it will be denoted as "contact laying" (or simply "laying"). However, the same regions may, instead, originate from another mechanism (undistinguishable from contact laying in Fig. 4): some tangential motion of nodes along the contact interface, due to local slip. In case the slip field would be heterogeneous, it could cause in-plane compression or dilation of the contact region, leading to variations of the contact area. Below, this third elementary mechanism potentially causing area variations (either increase or decrease) will be denoted as "in-plane deformation", under the form of either "in-plane compression" (or simply "compression") or "in-plane dilation" (or simply "dilation"). The existence of micro-slip in the model can be unambiguously identified in Fig. 6, where the evolution of the contact during incipient loading is shown at five select tangential displacements of the rigid plate, identified on Fig. 3(a): d 0 = 0, d s at the onset of gross sliding and three intermediate configurations (d i ). The nodes of the elastic surface are labeled in red if they have undergone some local slip, and in blue if they are still stuck to the same point of the rigid plate. It can be seen in Fig. 6 that the slip zone advances from the contact periphery towards the centre at the expense of the stick zone until the stick zone vanishes and full sliding occurs. The corresponding Movie S1 is provided as a supplementary material (shaded strips are introduced in the movie to visualize the in-plane deformation of the contact surface). The contact zone and the stick zone have an elliptical shape, the related shear-induced anisotropy increases with increasing displacement d and is higher for the lower normal loads, all effects also observed in the 8 Green regions indicate the current (deformed) location of the initial contact zone that has been lifted.
experiments . Partial slip configurations as those seen in Fig. 6 are classically found in models of sheared frictional linear elastic sphere/plane contact (Johnson, 1985). They typically correspond to heterogeneous slip fields, thus causing in-plane deformations. Figure 7 provides a detailed insight into this mechanism, by showing the field of local relative area changes. For each contact node, its tributary area both in the initial configuration (at d 0 ), A i 0 , and in the current configuration, A i , are computed. The color map in Fig. 7 corresponds to the relative area change, (A i − A i 0 )/A i 0 . In-plane dilation is observed close to the leading edge while compression is observed close to the trailing edge, which is consistent with the distribution of the σ xx stress component shown in Fig. A.16(c). The decreasing green elliptical region in the middle of the contact zone, in which the relative area change is equal to zero (green), corresponds to the stick zone in Fig. 6.
In order to quantify the respective contributions of the three possible mechanisms of area variations (lifting, laying and in-plane deformations), we used the following procedure. First, the contact nodes are grouped into six sets that are defined according to the contact state and relative contact area change. Set S 0 is determined at d 0 , the remaining sets are determined at each current displacement d. The following sets are defined: S 0 : nodes in contact at d 0 , i.e., after applying the normal load; S d : nodes currently in contact; S lift : lifted nodes, i.e., those in S 0 but not in S d ; S lay : new nodes in contact, i.e., those in S d but not in S 0 ; S comp : nodes in S 0 and in S d , for which the tributary area decreased, A i < A i 0 ; S dil : nodes in S 0 and in S d , for which the tributary area increased, A i > A i 0 . The total initial and current areas of the nodes belonging to each set are defined as a 0 (S) = i∈S A i 0 and a(S) = i∈S A i , so that, in particular, A 0 = a 0 (S 0 ) and A = a(S d ). Finally, the 9 contributions to the total relative area change that correspond to the various individual mechanisms are defined as:  Figure 8 shows how the individual mechanisms contribute to the total relative contact area change in the neo-Hookean model. For all normal loads, the model predicts that lifting represents about 90-95% of the total shear-induced area reduction, and is thus the primary mechanism responsible for it. The contributions of dilation and compression are individually significant, however, they essentially cancel out, so that in-plane deformation is only responsible only for about 5-10% of the total area reduction. Contact laying is rarely observed (if so, at the leading edge only), so that its contribution to the total area reduction is essentially negligible in the model.

Illustration experiment
In this section, we perform an illustration experiment, similar to those of Sahli et al. (2018Sahli et al. ( , 2019 used to validate the model in the previous section, in order to test whether the elementary mechanisms responsible for contact area reduction in the model (lifting, laying and in-plane deformation) can also be identified experimentally. The strategy is to incorporate particles close to the surface of the elastomer sphere and use them as tracers of the local motion of the frictional interface during incipient tangential loading. Figure 9: Sketch of the opto-mechanical setup. The lower glass plate supporting the elastomer sample (inset) is attached to an optical table, while the substrate (upper glass plate) is driven tangentially at constant velocity V, through a horizontal double cantilever. The tangential force along x is measured at the right extremity of the cantilever. The contact interface is illuminated from the top (lighting system not shown) and imaged in reflection with the camera. Inset: Sketch of the elastomer sample (blue, top). It has a cylinder-like shape, the top of which features a spherical cap. Particles are buried about 16 µm below the cap surface and serve as displacement tracers. The whole sample is attached to a glass plate (grey, bottom).  Fig. 9. The two main improvements consisted in (i) replacing the single beam cantilever with a double beam cantilever and (ii) adding force sensors to measure the normal load. This setup is used to shear the interface between a smooth glass plate and a cross-linked PDMS sample with a smooth, particle-seeded spherical cap (Appendix B.2). The experiment presented here was performed under constant normal force P = 1.85 N. 40 s after the contact has been created, a constant driving velocity V = 0.1 mm/s was imposed on the glass, over a total distance of 2 mm. The evolution of the tangential force Q as a function of the displacement of the glass plate is shown in black in Fig. 10. The incipient tangential loading of the interface reaches a maximum, denoted as the static friction force Q s , before a full sliding regime during which Q < Q s . Figure 10: Concurrent evolution of the tangential force Q (solid black curve, left axis) and of the contact area A (solid grey curve, right axis) as a function of the imposed tangential displacement of the glass substrate. P = 1.85 N. V = 0.1 mm/s. Four values of the displacement are indicated, which will be used in further figures: d 0 = 0 (blue, dotted), d s (red, dash-dotted) the displacement at the static friction peak Q s , and two intermediate displacements d 1 (cyan, dashed) and d 2 (yellow, solid).

Experimental methods
Typical raw images of the contact interface are shown in Fig. 11(a-d) for four different displacements: d 0 before any shear, d s the displacement at the static friction peak Q s , and two intermediate displacements d 1 and d 2 (shown in Fig. 10). A full movie is available as supplementary material (Movie S2). The contact region corresponds to the biggest region of dark pixels, because the light rays are transmitted through the PDMS/glass interface and absorbed by a black layer inserted at the bottom of the sample. Figures 11(e-h) show typical contours of the contact, the inner area of which defines the contact area, A. The evolution of A as a function of the displacement imposed to the glass plate is shown in Fig. 10 (grey curve). Note that both curves in Fig. 10 show the same qualitative behaviour as that in previous experiments in the literature (Waters and Guduru, 2010;Sahli et al., 2018;Mergel et al., 2019), indicating that the introduction of the particles into the elastomer sphere is not significantly affecting the mechanical response of the interface. The oscillations observable in the A(d) curve in Fig. 10 originate from so-called 12 re-attachment folds (Barquins, 1985) occurring at the trailing edge of the contact, and clearly visible in Movie S2. The lighter gray background corresponds to out-of-contact regions where a fraction of the light intensity is reflected back to the camera, at the interface between the bottom face of the glass substrate and the air. The small region of dark pixels seen on the top right of the contact is due to a black marking drawn on purpose on the top surface of the glass substrate. It serves as a tracer to monitor the macroscopic motion of the substrate. In (a-d), the main dark region is the contact, the bright spots are the particles, and the small dark region in the top right corner is a tracer drawn on the glass substrate. In (e-h), the contours are those of the contact (same colors as in Fig. 10), the inner area of which defines the contact area; the black spots correspond to the particles and are the objects used for the tracking procedure. Note that images (a-d) are in the frame of the camera, while (e-h) are in the frame of the (moving) glass plate.
The contact region is not uniformly dark, as it was in our previous sphere/plane experiments (Sahli et al., 2018Mergel et al., 2019), but is sprinkled with random bright spots which are due to the reflection of light on the particles that we introduced close to the elastomer surface (Appendix B.2). By using a homemade tracking procedure (Appendix B.3), we were able to use those spots as tracers of the contact evolution and to generate a dataset in which the position (T x (t), T y (t)) of each tracer, in the frame of the glass plate, is given for each image, i.e., for each time step. Note that only the tracers at the vertical of the contact regions are visible, so that they will be used not only to measure in-plane displacements but also to look for lifting (resp. laying), which would correspond to the disappearance (resp. appearance) of tracers at the periphery of the contact.

Results and analysis
In this section, based on the results of the tracking procedure, we perform a thorough analysis of the tracers' behaviour (appearance/disappearance and relative in-plane motion) to test the existence, in the experiments, of the elementary mechanisms found to be responsible for contact area evolution in the model (see Section 2.3). We first demonstrate that contact lifting and contact laying do occur, at opposite sides of the periphery of the contact zone, and we quantify their amount using Voronoi tessellation. Second, using Delaunay triangulation, we illustrate the progressive development of a heterogeneous in-plane strain field within the contact area.

Contact lifting
As already mentioned in Section 3.1, in the contact images, only the particles at the vertical of the contact region can be seen. This opens the possibility to check whether lifting occurs at the 13 interface. Indeed, the signature of local contact lifting is a particle that disappears from the image when reached by the moving contact periphery, meaning that a point of the elastomer initially in contact with the glass has been lifted out-of-contact. In practice, we looked for tracer trajectories that ended at a location closer than 20 pixels to the contact periphery. All such tracers are represented as disks in Fig. 12. Filled disks correspond to the location of the lifted tracers in the initial image (at d 0 ), open disks indicate the location of tracers when they disappear. It appears that a significant number of tracers are indeed lifted during the incipient tangential loading of the contact. The large majority of them are found at the trailing edge of the contact (at the right of Fig. 12, where points of the glass leave the contact). Two of them are found at the leading edge of the contact, but are probably tracers initially close to the contact periphery and disappearing due to noise in the images. Figure 12: Measurement of the lifted and laid areas. Solid colored lines: contours of the contact at the four selected displacements (same colors as in Figs. 10 and 11). Open (resp. filled) disks: position of the lifted particles when they disappear (resp. at d 0 ). The symbol color corresponds to the displacement at which the particle is lifted (same color code as for the contours). Grey cells: cells of the Voronoi tessellation at d 0 associated to lifted particles at d s . Filled (resp. open) diamonds: position of the laid particles at d s (resp. when they appear). The symbol color corresponds to the displacement at which the particle are laid. Red cells: cells of the Voronoi tessellation at d s associated to laid particles at d s .
The color of each disk in Fig. 12 corresponds to the glass displacement at which the tracer disappears. It appears that the colors of the filled disks are spatially organized, from blue (early disappearance) for the rightmost disks to red (late disappearance) for the leftmost disks. Such a pattern indicates that lifting occurs through an inward front propagation, starting from the trailing edge of the contact.
In order to quantify the amount of contact area that is lost due to lifting, we attempted to assign a representative elementary area to each tracer. For this, we performed, in the initial image, a bounded Voronoi tessellation on the centroids of the tracers. The boundary of this tessellation corresponds to the contact contour. Each tracer was thus assigned the area of its Voronoi cell. Then, at each subsequent image, the lifted area is computed as the sum of the areas of the Voronoi cells of all lifted tracers. The evolution of the lifted area along the experiment is shown in blue in Fig. 13. The evolution is stair-like, because each lifting event increases abruptly the lifted area by a finite amount (the cell area). The estimated final lifted area represents about 31% of the initial contact area. Figure 13: Contributions of the various mechanisms to the area variation. Ratio of area variation, ∆A, to initial contact area, A 0 , as a function of the imposed displacement, d. Black: reduction of the total contact area. Blue (resp. red): area lost (resp. gained) by lifting (resp. laying), measured as in Fig. 12. Green: area variations due to slip-induced in-plane deformation, measured as in Fig. 14.

Contact laying
The mechanism of laying, i.e., points of the elastomer which get into contact with the glass upon shearing, can be analyzed in a similar way as lifting. We performed the same analysis described in the previous section, but backward in time: the first image considered was actually the one at the static friction peak (at d s ), while the last one was the initial image (under zero shear, at d 0 ). Doing so, tracers appearing close to the contact periphery in forward time are detected as disappearing in backward time.
The results of this analysis are shown in Fig. 12. The symbols for laid tracers are diamonds, either filled for their position in the final configuration (at d s ) or open for their position when they appear in the image. All laid tracers are found at the leading edge of the contact. Their color, which evolves from blue to red, corresponds to the stage at which they appear and is again spatially organized, indicating that contact laying occurs through an outward front propagation starting at the leading edge of the contact. The amount of contact area gained via laying is estimated using a bounded Voronoi tessellation performed in the final contact area (at d s ). The laid area at a given glass displacement is counted as the sum of the areas of all laid tracers at that 15 instant. In Fig. 12, the red cells correspond to the area assigned to laying in the final image (at d s ). The evolution of the area gained thanks to laying as a function of the glass displacement is shown in Fig. 13 in red. Just like the lifting curve, it is stair-like, as cells appear abruptly with a finite area. The final amount of laid area is about 9% of the initial contact area. Unlike lifting, the laying mechanism seems to initiate only after a finite shear has been applied.

In-plane deformation
Interestingly, in the previous analyses of both contact lifting and laying, the filled and open markers have different positions with respect to the glass substrate. This is a clear indication that slip occurs between the two instants: lifting of a given elastomer point is preceded by slip, while laying of an elastomer point is followed by slip. Slip is found to be roughly parallel, but opposed to the glass motion, and occurs in both the leading and trailing regions of the contact.
The two above observations are consistent with a scenario of a micro-slip front propagating inward the contact region as shear is increased, which is classical in sheared sphere/plane contacts, either rough (Prevost et al., 2013) or smooth (Chateauminois et al., 2010). In this scenario, which is apparent in Movie S2 (supplementary material) and also observed in the model results of this paper (see Section 2.3), a peripheral slip region progressively invades the contact, replacing a shrinking central stuck region. From the combination of backward slip at the trailing edge and a stuck zone at the center of the contact, one expects the building-up of in-plane compression of the elastomer material in the trailing half of the contact. Symmetrically, in-plane dilation is expected in the leading half. The overall effect of both types of in-plane deformation must lead to a change in contact area that we attempted to estimate in the following way.
In the up-coming analysis, we only considered the tracers present in the initial image and that could be followed during the whole experiment, from d 0 to d s , thus excluding the lifted and laid ones. We performed a Delaunay triangulation to mesh all those tracers in the initial image, as shown as a gray network in Fig. 14. As shear increases, the tracers move relative to each other. Keeping the same mesh, we follow the relative changes of areas of all triangles of the initial Delaunay triangulation. The local in-plane deformation is thus estimated as this relative area change. The color of each cell corresponds to its relative area change with respect to the situation at d 0 (see colorbar). Colder (resp. warmer) colors mean in-plane compression (resp. dilation). Fig. 14 shows the evolution of the field of in-plane deformation as shear is increased. The color code (cold colors for compression, warm ones for dilation) provides a clear picture of those in-plane deformations. Before d 1 no deformation is observed, presumably because the slip front has not reached yet the fraction of initial contact covered by the Delaunay triangulation. Between d 1 and d 2 , a compressed region sets in at the trailing side of the contact; the absence of a detectable symmetrical dilated region on the leading side suggests that the propagation of the slip front is not axisymmetrical, as usually considered in models of the incipient shear loading of circular contacts (Barber, 2018). At d s , a heterogeneous in-plane deformation field is fully developed, with a large compression region on the trailing side and a smaller dilation region on the leading side, consistent with the effect of a slip zone having invaded the whole contact area. Those field measurement are also qualitatively consistent with previous measurements on a similar system, but made only along the central line of the contact in steady-state sliding (Barquins, 1985).
The evolution under shear of the sum of the areas of all triangles contained in the deforming Delaunay triangulation is shown in green in Fig. 13. This area remains essentially unchanged during the first half of the experiment, and then progressively decreases by up to about 7% of the initial contact area.

Uncertainties
The uncertainties of the above experimental measurements are mostly related to the finite surface density of useful tracers. We identified 263 tracers strictly larger than 4 pixels in the initial image, on which the tracking procedure has been applied. Out of them, only 31 were useful for the Voronoi analysis (i.e., corresponding either to lifted or laid areas) and 94 for the Delaunay analysis. The average Voronoi area of the useful cells was thus about 0.18 mm 2 , i.e., about 1.54 % of the initial contact area. In practice, due to the randomness of the tracers' locations, the individual areas varied from 0.2 % to 4.0 % of the initial contact area. Because the Voronoi cells are the smallest area units in our lifting analysis, the evolution rate of the estimated lifting area suffers from large fluctuations (see drops up to about 3-4 % of A 0 in the blue curve of Fig. 13). In addition, the final estimate of the lifted area is actually an overestimate because, for the latest lifted cells, the fraction of their area to the left of the tracer has not been really lifted yet. Both effects would be reduced with a larger density of tracers, yielding a smoother curve with a smaller amplitude. The very same discussion holds for our estimate of the laid area.
Concerning the estimate of in-plane deformation, one should note that the region probed by the Delaunay-based analysis is only a fraction of the target region (the part of the initial contact that never lifts nor lays during the experiment). This is apparent for instance on Fig. 14(a), with the white strip between the leading part of the contact contour and the Delaunay triangulation. One expects that an analogous strip is also lost on the trailing side, between the Delaunay triangulation and the lifted region. Those strips are due to the absence of tracers arbitrarily close to the contact contour. Because the largest strains are precisely found on the periphery of the Delaunay triangulation (see Fig. 14(d)), the estimate of the area lost by in-plane deformation may be subject to a significant error. The value of 7 % of area reduction due to this mechanism must then be taken with caution. Again, those uncertainties would be reduced with a larger surface density of tracers.

Discussion
The model presented in Section 2 has been shown to capture quantitatively the experimental results of Sahli et al. (2018Sahli et al. ( , 2019 about the anisotropic shear-induced area reduction of sphere/plane elastomer contacts for relatively large normal loads (Fig. 5). We emphasize that this agreement is obtained without any adjustable parameter, in the sense that virtually all model parameters are set by independent measurements on the same experimental system. Those parameters are: all geometrical features of the elastomer sample, the shear modulus of the elastomer treated as a hyperelastic neo-Hookean solid, and the contact shear strength of the PDMS/glass interface. Note that the elastic bulk modulus (in our nearly incompressible framework), the mesh size and the timescale in the regularized Tresca model are three purely numerical parameters, the precise values of which have been checked to negligibly affect our results. In our opinion, such an absence of adjustable parameter provides our model with a decisive comparative advantage with respect to competing models of shear-induced contact area reduction. In particular, the fracture-like adhesion-based models mentioned in the Introduction, in order to offer a good quantitative agreement with experiments Das and Chasiotis, 2020), require prescription of one or more mode-mixity functions, the shape and amplitude of which need to be finely tuned . Also note that our model naturally captures the anisotropic evolution of the contact shape, whereas most fracture-based adhesive models assume axisymmetry of the contact (see  for the only exception, to our knowledge).
The direct implication of our model is that the finite-deformation effects and the non-linear elasticity of elastomers are presumably the key ingredients explaining those experimental results, rather than viscoelasticity (as already argued qualitatively in Sahli et al. (2018)) or adhesion. Nevertheless, our model has only been applied to contact configurations submitted to relatively large normal load, in the newton range (Sahli et al., 2018), while other experimental studies have considered much smaller normal loads, in the millinewton range (Savkoor and Briggs, 1977;Waters and Guduru, 2010;Mergel et al., 2019). For those lighter loads, adhesive stresses may be of the order of, or even exceed, contact stresses, and non-linear elasticity may not be the dominant ingredient for shear-induced contact shrinking anymore. We suggest that identifying the normal load regimes in which non-linear elasticity or adhesion dominate the shear behavior of frictional sphere/plane contact is an important goal for future contact mechanics models. In the following, we will restrict our discussion to the high-load regime that we could explore in this study.
By comparing models assuming either small or finite strains, we have shown that non-linear elasticity is the single model ingredient responsible for the significant shear-induced contact area reduction observed in sphere/plane elastomer contacts (Fig. 3). Out of the many hyperelastic models suitable for elastomers, we have selected two models, and we have shown that the contact response may significantly depend on the model, the Mooney-Rivlin model predicting a visibly higher area reduction than the neo-Hookean model (Fig. 3), and thus also higher than observed in the experiment. Independent of the specific non-linear elastic model chosen, the finite-strain formulation naturally includes the exact geometrically non-linear kinematics. Even if the effects cannot be separated, it is possible that it is the geometric non-linearity that is actually responsible for the shear-induced contact area reduction.
The contact area behaviour of our main, neo-Hookean model has been fully validated using a competing finite-element model based on different regularization options, meshing and solver (Appendix A.6). This good agreement demonstrates the robustness of our model results to all implementation details. In addition, since both model approaches yield somewhat different shear stress distributions, we suggest that the contact area reduction phenomenon may be mainly sensitive to the average shear stress within the contact area, rather than to the precise distribution of the local shear stress.
Pinpointing the exact qualitative way in which finite deformations and non-linear elasticity amplify the contact area reduction still remains elusive. What we instead managed to do is to quantify, for the first time, the various elementary mechanisms by which area variations occurs. Only three of them can exist: contact lifting (elastomer initially in contact which is lifted out-of-contact), contact laying (creation of new contact) and in-plane deformation. Our model results suggest that, in the experimental conditions used in Sahli et al. (2018Sahli et al. ( , 2019, lifting is the dominant mechanism explaining the observed significant area reduction and its anisotropy (Fig. 8). This conclusion is found true for all normal loads explored, with the total reduction amplitude (and thus the lifting amount) being dependent on the normal load (as also observed experimentally). In a future work, it would be interesting to vary systematically all model parameters beyond the experimental range, in order to clarify their respective roles in the relative amplitudes of the three mechanisms potentially responsible for area variations.
Besides lifting, while laying is found negligible in the model for all normal loads, in-plane deformations are locally non-negligible. However, their overall contribution to area reduction remains modest because the area of the compressed regions decreases by about the same amount as the area of the dilated regions increases. Such in-plane deformations are fully related with the existence of a micro-slip front propagating inward from the contact periphery as the tangential load is increased (Fig. 6). Such a phenomenology is classical in linear elastic models of the incipient loading of frictional sphere-plane contacts (Johnson, 1985;Barber, 2018). In contrast, in our finite-strain elastic model, the micro-slip front clearly looses its circularity (see Fig. 6). Overall, the incipient shear-loading of smooth sphere/plane elastomer contacts is thus characterized in the model by two propagating fronts: a lifting front at the contact periphery, and a micro-slip front delimiting a shrinking-then-vanishing central stick region within the contact. The existence of those two different fronts, at different locations within the contact region, although explicitly acknowledged by some authors (Johnson, 1997;McMeeking et al., 2020), has never been properly described in fracture-based adhesive models. Indeed, the mode-mixity function used in those models is intended to describe the additional energy dissipation in the contact through frictional micro-slip, but such a dissipation is assumed to be located at the contact periphery, and not within a growing, finite region of the contact.
The relevance of all the above model observations have been qualitatively confirmed by the original experimental results of Section 3. Using an elastomer sphere seeded with tracers, we demonstrated that shear-induced anisotropic contact area reduction is indeed accompanied by all three possible mechanisms. Although a higher areal density of tracers would be desirable to reduce the measurement uncertainties, we could conclude that lifting is presumably the dominant area reduction mechanism, also in the experiments. Although experiments revealed a potentially larger contribution of the laying and in-plane mechanisms, a definitive quantitative comparison with the model is hindered by the finite experimental resolution. Also note that, due to the higher contact shear strength (0.53 rather than 0.41 MPa) and smaller Young's modulus (1.4 rather than 1.8 MPa) of the present experiments compared to those of Sahli et al. (2018Sahli et al. ( , 2019), our numerical model could not converge in those more severe conditions, thus impeding direct comparison between model and tracer-including experiments. Such a comparison is also an important challenge for future works.
Finally, the simplicity and generality of our model assumptions suggest that our results may 19 be also relevant to other systems than the elastomer sphere/plane contacts studied here. First, non-linear elasticity being a generic feature of soft materials, from gels to human skin, we expect it to be a likely mechanism for contact area reduction in all studies involving soft materials other than elastomers. It would thus be interesting to re-interpret recent experiments like those of e.g. Das and Chasiotis (2020) on polyacrylonitrile or those of Sahli et al. (2018); Sirin et al. (2019) on human fingertips, from the standpoint of finite-deformation mechanics and non-linear elasticity. Second, it has been argued in Sahli et al. (2018Sahli et al. ( , 2019 that the mechanisms of shearinduced anisotropic contact area reduction may be the same in millimetric sphere/plane contacts and in individual micrometric junctions within rough contact interfaces. Such a similarity across scales suggests that the conclusions of the present work may also be used to further understand the shear behaviour of soft material multicontacts. In particular, one may expect non-linear elasticity to be a necessary ingredient when modelling such systems. This may in particular be an explanation for the fact that, in Scheibert et al. (2020), a multi-asperity model based on linear elasticity failed to reproduce quantitatively the multicontact results of Sahli et al. (2018).

Conclusion
We have developed the first non-adhesive, non-viscous model of shear-induced contact area reduction in soft materials. Quantitative agreement with the recent experimental results of Sahli et al. (2018Sahli et al. ( , 2019 on sphere/plane elastomer contacts has been obtained, without adjustable parameter, using the Tresca friction law and the neo-Hookean hyperelastic model. We have demonstrated that the necessary ingredients for this agreement are finite deformations and nonlinear elasticity.
A detailed analysis of the model has revealed the elementary mechanisms responsible for contact area reduction. Local contact lifting is the dominant mechanism, while in-plane dilation and compression, although significant, approximately cancel each other. Those predicted mechanisms, and their relative contributions to area reduction, have been confirmed experimentally by tracking tracers introduced in an elastomer sphere in contact with a glass plate, and applying an original analysis to the tracking data. The challenges associated with both the modelling and experimental approaches have been thoroughly discussed.
All our results suggest a new perspective on the phenomenon of shear-induced contact area reduction in soft materials. This currently highly debated topic has been dominated by interpretations based on a dominant role of adhesion. Here, instead, we suggest that finite-deformation effects and non-linear elasticity can be equally important, all the more so as large normal loads are considered. Clarification of the respective validity domains of finite deformations and adhesion remains as a major open issue on the topic. rations of the body are considered: the reference configuration Ω and the current (deformed) configuration ω. The deformation that brings Ω to ω is described by the deformation mapping ϕ such that x = ϕ(X), where X ∈ Ω and x ∈ ω denote the position of a material point in the respective configuration. Further, the deformation gradient F = Grad ϕ is defined, where Grad is the gradient with respect to the reference configuration Ω.
Equilibrium equation in the strong form is written in the reference configuration Ω, and, in the absence of body forces, reads where P is the first Piola-Kirchhoff stress tensor, Div is the divergence operator relative to the reference configuration Ω, and the elastic strain energy W = W(F) specifies the constitutive response of a hyperelastic body. The reference material model used in this work is the nearly-incompressible isotropic neo-Hookean model specified by W = W nH , where W vol (J) is the volumetric part of the elastic strain energy, J = det F, andĪ 1 = trb is the invariant ofb = J −2/3 b, b = FF T . Material properties are here specified by the shear modulus µ and bulk modulus κ.

Appendix A.2. Contact with a rigid plane
The contact formulation briefly described below is limited to the case of contact with a rigid plane, which is sufficient for the purpose of this work, see e.g. Wriggers (2006) for a more general presentation. The orientation of the plane is specified by the outward unit normal ν, and it is assumed that the motion of the rigid plane is restricted to a translation, thusν = 0.
Contact kinematics is specified by the normal gap g N and tangential slip velocity v T that are defined for each point x on the potential contact surface γ c = ϕ(Γ c ), where Γ c denotes the contact surface in the reference configuration. In the present case, the kinematic quantities are simply given by where x R denotes the current position of a fixed point on the rigid plane, ⊗ denotes the diadic product, and (1 − ν ⊗ ν) is the projection operator on the tangent plane. The contact traction vector t is defined as the traction vector exerted by the body on the rigid plane. By the action-reaction principle it is equal to the traction vector acting on the body with a minus sign, and thus we have t = −σn, where σ is the Cauchy stress tensor, n is the outward unit normal to the contact surface γ c , and whenever contact occurs we have n = −ν. Note that t is a spatial traction vector, i.e., it refers to a unit area in the current configuration. The contact traction is decomposed into its normal and tangential parts relative to the rigid-plane normal ν, namely With the definitions introduced, the unilateral contact conditions can be written in the following standard form, The friction model is discussed in Appendix A.3.
To conclude the contact formulation, the virtual work principle, i.e., the weak form of the mechanical equilibrium, which is the basis for the finite-element implementation, is provided for the problem at hand, viz.
where δϕ is the test function that vanishes on the part of the boundary of Ω on which the displacement is prescribed, and we have δg N = ν · δϕ and δ g T = (1 − ν ⊗ ν)δϕ (cf. Eq. (A.4)).
Since t N and t T are spatial tractions, the contact contribution, i.e., the second term in Eq. (A.7), is integrated over the contact surface γ c in the current configuration.

Appendix A.3. Regularized Tresca friction model
In the Tresca friction model, the limit friction stress, denoted here by σ and called the contact shear strength, is constant and independent of the normal contact traction t N (Fig. A.15(b)). The Tresca model can be written in the following form: and the equations above hold in the case of active contact, i.e., for g N = 0. In the case of separation, i.e., for g N > 0, the tangential traction vanishes, t T = 0. Note that, in view of the unilateral contact condition (Eq. (A.6)), penetration (g N < 0) is not allowed. The first condition in Eq. (A.8) is the limit friction condition that constrains the norm of the tangential traction t T not to exceed the contact shear strength σ. The second condition is the slip rule that implies that the direction of the tangential (slip) velocity v T is that of the tangential 22 traction t T . Finally, the third (complementarity) condition controls the stick/slip state such that sticking contact (v T = 0) occurs when t T < σ and sliding ( v T > 0) may occur only when t T = σ.
One of the features of the Tresca model, as presented above, is that the tangential traction t T may suffer a jump at the instant of transition between active contact and separation. This, in turn, may lead to significant convergence problems in the respective computational scheme, for instance, in the framework of the finite-element method. Accordingly, a regularization scheme has been developed, as described below, and employed in the actual computations. Note that an alternative approach, based on the Coulomb-Orowan friction model (Fig. A.15(c)), has also been examined, and the respective results are briefly presented in Appendix A.6; the Tresca model can be interpreted as a limit case of the Coulomb-Orowan model for f → +∞.
In the regularization scheme adopted here, it is assumed that the tangential traction t T does not drop to zero instantly after separation occurs but requires some characteristic time τ to gradually vanish. Accordingly, the following simple evolution law is adopted in the case of separation so that the magnitude of the tangential traction t T decreases towards zero at a constant rate of σ/τ. A similar regularization, leading to a kind of rate-and-state friction law, has been adopted by Cochard and Rice (2000) to regularize the abrupt changes in friction associated with abrupt changes in the normal contact traction, see also Prakash and Clifton (1993); Prakash (1998) for the experimental justification and respective constitutive modelling. Here, the regularization is introduced purely for computational reasons.

Appendix A.4. Finite-element implementation
For the finite-element discretization of the nearly-incompressible hyperelastic solid at hand, the eight-node hexahedral TSCG12 element (Korelc et al., 2010) is used. The element employs the enhanced assumed strain (EAS) formulation to circumvent volumetric and shear locking effects and has proven to perform very well in the contact problems considered in this work. As a means of verification, selected simulations have been repeated using the F-bar element (de Souza Neto et al., 1996), and the effect of the finite-element formulation has been found negligible.
Consistent with the discretization of the solid, the contact surface is discretized into four-node quadrilateral segments. Nodal quadrature is applied to the contact contribution in the virtual work principle (Eq. (A.7)) so that, effectively, the contact and friction conditions are evaluated at the individual nodes of the contact surface. In this formulation, the tributary area is determined for each contact node (cf. Lengiewicz et al., 2011), and the changes in the contact area can be traced by summing up the tributary areas of the nodes that are in contact at a given instant. Recall that the contact contribution is evaluated in the current configuration, and so is the nodal tributary area.
The unilateral contact and friction conditions (Eqs. (A.6) and (A.8)) are enforced using the augmented Lagrangian method (Alart and Curnier, 1991). In this method, the contact Lagrange multipliers are introduced as global unknowns, and the resulting system of non-linear equations is solved simultaneously with respect to the nodal displacements and Lagrange multipliers using the semi-smooth Newton method. As a result, the contact conditions are enforced in a numerically exact manner.
The beneficial feature of the augmented Lagrangian method is that the contact/separation and stick/slip states are uniquely defined at each global Newton iteration when the contact conditions (Eqs. (A.6) and (A.8)) are not yet satisfied. In the time-discretized setting, a special treatment is, however, necessary when switching between the Tresca friction model (Eq. (A.8)) and its regularization (Eq. (A.9)). In the direct approach resulting from the implicit time-integration scheme, the current contact/separation state determined at each iteration is used to choose whether the Tresca model (in the case of contact) or the regularization scheme (in the case of separation) is employed. However, we have observed that this leads to severe convergence problems. Accordingly, solely with the aim to determine whether to use the Tresca model or its regularization, the contact/separation state is here taken from the last converged time step. Furthermore, an adaptive time-stepping procedure is applied, which conveniently alleviates the difficulties caused by the high non-linearity of the problem.
Computer implementation has been performed using AceGen, a code generation system that employs the automatic differentation (AD) technique (Korelc, 2009;Korelc and Wriggers, 2016), see also Lengiewicz et al. (2011) for the details of the automation of the finite-element contact formulations. The computations have been performed in AceFEM, a flexible finite-element environment that is interfaced with AceGen.

Appendix A.5. Finite-element model
The geometric parameters and material properties used in the simulations are directly taken from one of the experimental datasets reported in Sahli et al. (2018) (smooth PDMS-sphere/bareglass contact, see Fig. 2C therein). The boundary conditions and the loading program also correspond to those in the experiment. The displacements are fully constrained at the bottom surface of the elastomer sample (cf. Figs. 1 and 9), except for a rigid-body translation in the z-direction so that the normal load P can be applied against a rigid plane representing the upper glass plate. Subsequently, a constant tangential velocity of the rigid plane is prescribed along the x-direction while the normal load P is kept constant. Further, the symmetry with respect to the y = 0 plane is exploited so that only one half of the sample is effectively modelled with adequate boundary condition imposed on the symmetry plane.
The resulting finite-element mesh is shown in Fig. 2. The sample is discretized with about 100,000 hexahedral elements, which gives approximately 380,000 degrees of freedom. The hanging-node technique is used to conveniently refine the mesh in the vicinity of the potential contact surface, and the size of the refined-mesh region is adjusted to the normal load. The contact surface comprises about 12,600 nodes with the size of the quadrilateral elements ranging from 13 µm for P = 0.27 N to 25 µm for P = 2.12 N.
The elastic properties of the elastomer sample have been identified by matching the contact area obtained from the computational model with neo-Hookean elasticity for the highest normal load (P = 2.12 N) and zero tangential load (Q = 0) to that measured in the experiment. Assuming a nearly incompressible response with the Poisson's ratio ν = 0.49, the Young's modulus has been identified as E = 1.80 MPa, which corresponds to a shear modulus µ = 0.60 MPa. The relationship between the parameters (E, ν) and (µ, κ) (cf. Eq. (A.2)) is given by µ = E/(2(1 + ν)) and κ = E/ (3(1 − 2ν)). It can be seen in Fig. 5 that, for the shear modulus determined as described above, the contact area at zero tangential load is correctly represented in the whole range of normal loads for both the neo-Hookean and Mooney-Rivlin models.
The value of the contact shear strength σ = 0.41 MPa in the Tresca model (Eq. (A.8)) is taken as the slope of the best linear fit of the evolution of the static friction force, Q s as a function of the concurrent contact area, A s = A(Q = Q s ) in the experimental dataset used here as a reference.
The Prakash-Clifton-like regularization parameter in the Tresca model is assumed to be τ = 0.2 s, a value which has been checked to be small enough to have negligible influence on the results. The computations have been carried out using an adaptive time stepping procedure. The typical number of resulting time increments was between 50 and 100, which corresponds to an average time increment of the order of 0.1 s. Note that the data points included in the figures reported throughout Section 2 and Appendix A correspond to selected instants. Sample results of finite-element computations are provided in Fig. A.16. The maps of three selected components of the Cauchy stress tensor are shown in the deformed configuration corresponding to full sliding at P = 2.12 N. The σ zz component is shown in Fig. A.16(a). Its value at the contact surface corresponds to the normal contact traction, and the corresponding distribution has a typical Hertz-like appearance with the lowest value in the middle of the contact zone (recall that the normal contact traction t N is negative). The shear stress σ xz is shown in Fig. A.16(b) with a constant value at the contact surface, in agreement with the Tresca friction law. Finally, the σ xx stress is shown in Fig. A.16(c) featuring zones of tensile stresses next to the leading edge and compressive stresses next to the trailing edge. The zones of tensile and compressive stresses correspond to the zones of surface dilation and compression, respectively, as discussed in detail in Section 2.3.
Note that the contact problem under consideration is highly demanding from the computational point of view. In the spatially discretized setting, the Tresca model may lead to abrupt changes of the contact nodal forces, and the regularization scheme (Eq. (A.8)) has been introduced to overcome the related convergence problems. Moreover, the frictional shear strength σ is of the same order as the elastic shear modulus µ, and the Tresca model implies that the corresponding high tangential tractions are applied step-wise at the contact zone boundary. This may lead to a significant distortion of the finite-element mesh, and the boundary of the contact zone exhibits the related mesh-dependent features, as illustrated in Fig. A.17(a). Mesh-dependent features in the form of secondary contact regions are also observed, see Fig. A.17(b), particularly at lower normal load P and higher tangential load Q (cf. Fig. 6). Since those separated contact spots have a small total area compared to that of the main contact region, they have been included in all evaluations of the total contact area, A. In contrast, when we evaluated the contact sizes along 25 and orthogonal to shear, and ⊥ respectively, those separated spots have been excluded because they would have induced a significant noise and error, especially on (see e.g. Fig. A.17(b)). Due to the high computational cost and significant convergence challenges, we were not able to solve the problem for a finer mesh. However, we have carried out additional simulations for a twice coarser mesh, and we have checked that the results are essentially identical to to those obtained for the reference mesh and thus are not affected by the mesh-dependent features discussed earlier. This, in particular, concerns the evolution of the contact area as a function of the tangential force, as shown in Fig. A.18.
Appendix A.6. Alternative model using Coulomb-Orowan friction To test the robustness of the results obtained using the time-regularized Tresca friction model (Eq. (A.9)) additional computations have been carried out using the Coulomb-Orowan model (Fig. A.15(c)) as an alternative regularization approach. The friction coefficient in the Coulomb-Orowan model has been adopted equal to f = 10, which in practice was the highest achievable value leading to convergence. All geometrical and material parameters and boundary conditions are otherwise the same as in the main (Tresca-based) model.
The unilateral contact and friction conditions were enforced using the penalty method. In particular, the tangential tractions could reach the contact shear strength, σ only when a slip distance of 0.3 % of the contact size was reached. Such an elastic slip allowed to smoothly reach the final shear state and to converge faster. Using this approach, the evolution of the contact area versus the tangential load can be consistently extracted, where the intermediate states correspond to a progressive transition of the contact zone from the sticking to sliding state.
The corresponding computations have been carried out using ABAQUS 2020 finite-element package. The hyperelastic (neo-Hookean) solid was meshed with 8-node linear brick hybrid elements with reduced integration (C3D8RH). In order to alleviate the issues related to the almost incompressible nature of the material and volumetric locking problems, the hybrid formulation was used (Nguyen et al., 2011). The mesh was refined down to 30 µm in the vicinity of the contact zone for all normal loads, leading to a total of nearly 145,000 hexahedral elements and about 580,000 degrees of freedom, see Fig. A.19(a).
The large-displacement formulation was employed to account for the various non-linearities and solved using an implicit integration scheme based on the Newton-Raphson method. To overcome convergence problems, the simulation was perfomed in two steps. In the first step, the normal load P was applied in frictionless conditions. In the second step, the sliding velocity and the friction conditions were applied.
From a qualitative point of view, the predicted shape of the contact morphology shown in Figure A.19(b) for the full sliding configuration, is found consistent with that of the Tresca model (compare to Fig. 4(b)). An anisotropic change was also observed with a stronger reduction at the trailing edge of the contact. Figure A.19(c) summarizes the quantitative evolution of the contact area for the same four loading cases considered in the Tresca model. The predicted area reduction is in good agreement with the experimental data and thus also with that obtained using the Tresca model ( Fig. 5(a)). Such a good correlation shows that the main model results are robust against the details of the regularization method adopted.
The main limitation of the Colomb-Orowan model is that the contact shear strength σ cannot be reached at the very periphery of the contact, where the local pressure remains smaller than σ/ f . In this region, the limiting shear stress is set by the friction coefficient f rather than by the contact shear strength σ. As a consequence, the maximum tangential force accessible in the simulations remains smaller than the product of σ with the contact area, and thus the simulated curves do not reach perfectly the desired maximum tangential force, denoted as the red line in tine was used to record the acquisitions of the forces, P and Q, and MATLAB (The MathWorks, Inc.) software was used to synchronize the data with the images.

Appendix B.2. Elastomer and glass samples
The particle-filled cross-linked PolyDiMethylSiloxane (PDMS, Sylgard 184, Dow Corning) sample with a spherical cap was prepared as follows. The elastomer base and curing agent are respectively mixed in a 10:1 weight ratio and degassed in a vacuum chamber to remove air bubbles introduced during mixing. Cross-linking is performed in two steps. First, a drop of PDMS mixture is poured into a plano-concave glass lens with a radius of curvature of 9.42 mm (Edmund Optics), spin-coated (SPIN 150, Spincoating) for 30 s at 2500 rpm and cured in an oven (Prolabo, Astel S.A.) at 80 • C for 20 min, to obtain a thin cross-linked PDMS layer, approximately 16 µm thick. The exposed surface of this layer is then sprinkled with silver particles, obtained from drying an Electrodag 1415M (Agar Scientific Ltd) diluted in acetone. In a second step, the lens with the particle-seeded PDMS layer is placed at the bottom of a cylindrical cavity (30 mm in diameter), which is then filled with the PDMS mixture, covered with a glass plate (40 mm×35 mm×5 mm) and cured at room temperature for 48 h. Once demolded, the sample has a 6 mm-high cylinder-like shape, the top of which features a spherical cap of a radius of curvature 9.42 mm, with a summit 2 mm higher than the top surface of the cylinder (see inset of Fig. 9).
A 5 mm-thick smooth bare-glass plate (Mirit Glass) was used as a slider. Before experiments, the glass surface was first gently pre-cleaned mechanically with distilled water using a lens cleaning tissue. It was then left in three different ultrasonic baths, each for 15 min, in the following order: soapy water (Decon90), ethanol and distilled water. Between the baths, the plate was rinsed in distilled water. It was finally dried in an oven at 90 • C for 15 min. Appendix B.3. Image analysis: segmentation and particle tracking The first step in the image analysis is the segmentation of the raw images into in-and out-ofcontact pixels. For this, we first crop the images around the contact, thus excluding the marker on the glass substrate (see top right black spot in Figs. 11(a-d)). Then, we transform the raw grayscale cropped images into binary images by thresholding. The threshold value is kept constant for all images and is selected by applying Otsu's method (Otsu, 1979) to the initial image. The segmented images can be directly used to measure the contact area, A, through selection of the largest dark object in the image and filling up of its inner holes, due to the presence of the particles. Typical contours of the contact are shown in .
The second step consists in the following denoising procedure applied to all images: (i) all objects whose area is a single pixel are removed, then (ii) the closest remaining objects are connected by successively performing one dilation and one erosion operation with a 3×3 square structuring element.
During the third step, we identify, in the initial image, all the objects whose area is strictly larger than 4 pixels. Those objects from the initial image are now considered as the tracers of the contact evolution for upcoming analysis. Note that for each tracer, we store the x-and y-widths of the smallest rectangle enclosing the tracer (W x and W y ).
The fourth step is to find the trajectory of each tracer by searching its successive positions in all images. For each couple of two consecutive images, the following tracking procedure is applied for each tracer: • extract two sub-images (first and second) from both consecutive images: their common position corresponds to the location of the tracer in the first image and their common widths correspond to W x and W y increased by 5 pixels on all four sides, to conservatively account for the dilation operation that will be performed in the next step. This 5 pixels widening is way larger than the distance (∆d max ) traveled by the slider between two images in steady sliding regime, ∆d max = 1 µm ≈ 0.14 pixel.
• perform a dilation operation with a 3×3 square structuring element on the first sub-image.
• multiply pixel by pixel the two resulting sub-images, i.e., the dilated first one with the unmodified second one.
• find all objects (connected components) in the multiplication image. Store the position of the object as the coordinates (T x , T y ) of its center of mass in the full image. In the case of multiple objects, only the closest object to the center of the sub-image is considered. If there is no object, then the trajectory ends.
The fifth step consists in removing all the trajectories that contain a non-realistic change, greater than 3 pixels, in the position of a tracer between two successive images.
In the last step, we first measure the glass plate displacement by tracking the position of the centroid of the tracer drawn on the glass for each image. We then translate the segmented and denoised counterparts of the raw images ( Fig. 11(a-d)) in such a way that this tracer remains at a constant position. In other words, we place our images in the frame of the glass substrate ( Fig. 11(e-h)), which actually moves from left to right at velocity V in the laboratory frame. Finally, we subtract the glass plate displacement to all our trajectories.
The outcome of the tracking analysis is thus a dataset in which the position (T x (t), T y (t)) of each tracer is given in the glass frame, for each image (i.e., for each time step).