Analytical and numerical results for the elasticity and adhesion of elastic films with arbitrary Poisson’s ratio and confinement

ABSTRACT We present an approximate, analytical treatment for the linearly elastic response of a film with arbitrary Poisson's ratio , which is indented by a flat cylindrical punch while resting on a rigid foundation. Our approach is based on a simple scaling argument allowing the vast changes of the elastomer’s effective modulus with the ratio of film height and indenter radius to be described with a compact, analytical expression. This yields exact asymptotics for large and small reduced film heights , whereby it also reproduces the observation that has a pronounced minimum for at . Using Green’s function molecular dynamics (GFMD), we demonstrate that the predictions for are reasonably correct and generate accurate reference data for effective modulus and pull-off force. GFMD also reveals that the nature of surface instabilities occurring during stable crack growth as well as the crack initiation itself depend sensitively on the way how continuum mechanics is terminated at small scales, that is, on parameters beyond the two dimensionless numbers and defining the continuum problem.


Introduction
This paper revisits the contact mechanics of confined, linearly elastic layers of height h sandwiched between a rigid surface and a circular rigid punch of radius a. A central quantity of such films is the effective modulus � E [1] as a function of the reduced height h=a and the Poisson's ratio ν. � E is defined as the ratio of mean contact stress and relative height change. The arguably most important reason for wanting to know � Eðh=aÞ is that it allows the pull-off stress σ p � [2,3] and the fracture mechanisms of confined elastomers to be determined. [4][5][6][7][8] For large h=a or small ν, adhesive contact failure is sudden, i. e., the tensile force drops discontinuously from its maximum value to zero under quasi-static loading, even when the system is displacement driven. However, stable crack growth occurs for (nearly) incompressible elastomers once h=a falls below approximately two. As h=a drops below unity, the contact shape is no longer circular during stable crack growth but clearly symmetry broken [4,[6][7][8][9] so that contact features have a characteristic linear dimension λ minimizing the total energy. [4,8,[10][11][12][13] The limits of unconfined (h=a ! 1) [14] and incompressible (ν ¼ 0:5), highly confined (h=a ! 0) [15] elastomers were solved more than half a century ago. Over time, amendments to the latter case were made with respect to boundary conditions, the geometry of the elastomer, and other details. [1,4,12,16] However, intermediate confinements h=a � 1 have not yet been solved analytically for extended two-dimensional films, although solutions for poker-chip specimens [17,18] as well as elastomeric strips [19] can be found. Extended films have only been treated numerically with finite-element (FE) simulations, [2,3,5,7,20,21] which lead to the suggestion of semi-empirical relationships between � E=E � and h=a, [3,5,20] where E � ¼ E=ð1 À ν 2 Þ is the contact modulus and E the Young's modulus. They turn out to benefit the interpretation of real-laboratory experiments, in particular, to explain different crack propagation mechanisms during detachment. [4,[6][7][8] A deeper understanding of the � Eðh=aÞ dependence might also prove useful in interpreting observations made on confined elastomers in contact with rough indenters. [22][23][24] Unfortunately, most existing semi-empirical � Eðh=aÞ relations were only designed for Poisson's ratios equal to or just below 0.5 so that confinement effects of various soft materials with small Poisson's ratios, such as foams, corks or soft isotropic metamaterials with negative Poisson's ratio are not quantitatively understood. Moreover, as demonstrated in this work, the extreme confinement limit, in which deviations from ideal incompressible matter have not yet been described satisfactorily for the given elastic film geometry. While instabilities of the elastomer surface were studied analytically for arbitrary Poisson's ratios, [10,25] the focus was quickly laid on ideally incompressible elastomers. [10,12,21] In addition, no physically motivated, closed-form expressions for the dependence of the effective modulus on the reduced film height have been proposed.
The original main motivation for this article was to identify a nonempirical relationship for � Eðh=a; νÞ, which allows us to easily rationalize the minimum in � Eðh=a; νÞ and estimate the range of Poisson's ratios, in which an elastomer film behaves as if it were incompressible. To this end, we propose that the energy needed to deform the elastic film should be most sensitive to the stiffness of a surface undulation with a wavevector in the order of the inverse punch radius. Even if such a simple scale argument may not outperform existing, more empirical models for all possible combinations of h=a and ν, it should improve our ability to estimate � Eðh=a; νÞ, in particular in the limit of extreme confinement and/or small Poisson's ratio. It certainly behoves us to examine numerically the accuracy of any scaling relation, which we do by running Green's function molecular dynamics (GFMD) simulations. [26] This also allows us to produce reference data for the pull-off stress as a function of h=a and ν. While simulating the detachment process, we realized that the analysis of surface instabilities that occur during stable crack growth at h=a , < 1 is interesting in its own right. We therefore include an in-depth analysis of how substrate symmetry, lattice trapping, and stochastic irregularities in the form of thermal noise, as well as their interplay affect the patterns that occur when the surface morphology becomes unstable during detachment.
The remainder of this article is organized as follows: Model and numerical methods are introduced in Section 2. Our scaling approach is presented in Section 3. Section 4 contains a comparison between theory and simulations as well as additional simulation results. Conclusions are drawn in Section 5.

Model
The investigated model system consists of isotropic, linearly elastic films of varying film height h resting on a perfectly flat and perfectly rigid foundation with a surface normal in the z direction. The in-plane extent of film and foundation are taken to be infinitely large and a no-slip condition is assumed between them. The opposite surface of the elastomer interacts with a rigid circular punch of radius a through a hard-wall constraint with a slip condition. Such systems can be effectively simulated by assuming periodic boundary conditions in the xy plane, as long as the linear dimension L of the periodically repeated simulation cell exceeds the punch radius a by a sufficiently large padding, which is most effectively chosen to be larger than but of order minðh; aÞ.
In the just-defined setup, the elastic energy of the elastomer is given by [27][28][29] where q is an in-plane wave vector with absolute value q, ũðqÞ is the Fourier coefficient of the displacement field of the elastomer's surface facing the indenter and For an infinitely large system without periodic boundaries, the sum on the r.h. s. of Eq. (1) will be replaced with an appropriate integral. For some of the calculations presented in this study, knowledge of the asymptotes of cðν; qhÞ is useful. A Taylor expansion reveals them to be cðν; qhÞ ¼ 1 for qh � 1 c 1 ðνÞ=ðqhÞ for qh � 0:5 À ν and ν < 0:5 1:5=ðqhÞ 3 for qh � 1 and ν ¼ 0:5: with In our analytical treatment, the interaction between indenter and elastomer is a non-overlap constraint. In addition, a surface energy γ is gained per unit area where surfaces touch. The model is then replaced with a cohesive zone model for the numerical solution of the contact problem, which is described next.

Methods
The contact problems were solved numerically using Green's function molecular dynamics (GFMD) simulations. [26] GFMD is a boundary-value method, in which Newton's equations of motion for the displacement fields are solved in their Fourier representation. In compression simulations, we use an exact non-overlap constraint in conjunction with the fast-inertial relaxation (FIRE) [30] algorithm as described in Ref. [31] Typical simulations assume the linear dimension of the periodically repeated simulation cell to be three times the punch diameter and a discretization of the displacement field into 2; 048 � 2; 048 elements. While exploiting the circular symmetry of the problem would have allowed us to reduce the computational cost of the simulations substantially, we found it more time effective to use the implemented methods.
Although knowledge of V ela as a function of normal displacement and h=a determined from purely repulsive experiments is sufficient to deduce the adhesive pull-off force, see Section 3.3, simulations mimicking tensile loading were also conducted. This was not only done to double-check our pull-off force calculations but also to investigate the dynamics and failure mechanisms that occur during the detachment of confined elastomers. For this purpose, adhesion is modeled with a cohesive zone model (CZM), in which the gapdependent surface energy has the form [32] γðgÞ ¼ À γ � fcosðkgÞ þ 1g=2 for 0 � kg � π=Δa and zero else, where g is the gap between elastomer and punch. The parameter k was generally set such that the maximum stiffness of surface undulations was slightly more than twice the maximum (negative) curvature of the potential defining the CZM, i.e., k 2 ¼ 0:4 r q max E � cðν; q max hÞ, where q max ¼ ffi ffi ffi 2 p π=Δa and Δa the linear mesh discretization. In this way, the interaction is effectively as short ranged as possible while avoiding lattice trapping. The latter refers to a situation, where an individual degree of freedom, e.g., a GFMD discretization point, can have two or more mechanically stable positions, while all other points remain fixed. When addressing lattice trapping, the parameter k was set to 3.75 times its default value. To improve the convergence rate, the massweighting GFMD variant was used for adhesive simulations. [31] Computing time is furthermore reduced by progressively increasing spatial resolution and decreasing the rate of retraction upon approaching the point of maximum tensile force.
To also model the response of elastomers to small perturbations, some simulations were conducted at finite temperature with the help of a recently introduced GFMD thermostat. [33] To this end, the thermal energy was kept at about 0.1% of the adhesive energy gained in a single mesh element in which the elastomer makes perfect contact with the indenter.

Finite-size corrections
If the origin of the coordinate system coincides with the center of the flat punch, the macroscopic displacement u 0 is defined as uð0Þ À uðr ! 1Þ. Thus, the best simple estimate for u 0 when using a finite square-shaped simulation cell with length L is to replace u 1 ;uðr ! 1Þ with uðL=2; L=2Þ.
To reduce the finite-size error, we use a correction appropriate for semiinfinite elastomers [34] but damp it with the weight function wðh=aÞ ¼ tanhðh=aÞ for small h=a: u 1 � uðL=2; L=2Þ þ 5 wðh=aÞ fuðL=2; L=2Þ À uðL=2; 0Þg: The damping of the usual correction is needed because for finite h=a, the displacement field approaches u 1 exponentially quickly with increasing distance r from the origin rather than with 1=r.

Theory
As has been done before, [2,3,35] we define the effective modulus � E as the ratio of the mean (compressive) contact stress � σ ¼ F=ðπ a 2 Þ and the relative height change of the elastomer to the contact area, � ε ¼ u 0 =h, even if the uncompressed elastomer is a film rather than a free-standing cylinder of radius a, for which � E was originally introduced. [35] Here, u 0 is the normal displacement of the elastomer's surface right below the punch from its equilibrium height in the absence of an indenter. Thus, � E is given by Since a does not change with u 0 for a flat punch under compression, it follows that F is proportional to u 0 within linear elasticity. Thus, the elastic energy is simply given by As already argued in the introduction, the only in-plane length defining the contact problem is the punch radius a. Thus, under compression, dimensional analysis suggests that the elastic energy should predominantly reside in undulations with wave numbers of the order of q a ¼ 2π=a given that q ! 0 or q ! 1 modes are not dominant. In this case, a good estimate for the elastic energy would be In order to eliminate the big-O notation in Eq. (9), we introduce two proportionality factors α and β. Comparing the resulting elastic energy to Eq. (8) yields which is the central analytical result of this work. Sections 3.1 and 3.2 are concerned with a parametrization of α and β, while Section 3.3 summarizes how to deduce depinning forces from the h=a dependence of � E.
To further motivate our approach, Figure 1 shows displacement and stress fields in real space for the various confinements, which range from small (h=a � 1) and intermediate (h=a ¼ 1) via large (1 À 2ν � ðh=aÞ 2 � 1) to extreme (ðh=aÞ 2 � 1). These data are complemented by the elastic energy associated with individual ũðqÞ modes in the right column of Figure 1, which clearly supports the scaling hypothesis that energy predominantly resides in modes with wavelengths of order a.
The top row of Figure 1 shows the well-established properties of the flat-punch solution for semi-infinite elastomers. On the other end, in the two bottom rows, Figure 1 reveals a qualitative difference between large and extreme confinement, which may often be underappreciated, although the principle is known from works addressing poker-chip and elastic-strip geometries. [17][18][19] Specifically, for 1 À 2 ν < ðh=aÞ 2 � 1, the stress profile at the origin is close to an inverted parabola while for ðh=aÞ 2 � 1 À 2 ν, the stress is constant within (most of) the contact. It should also be noted that in all cases, a stress singularity occurs at the contact edge, which in our numerical treatment -and in reality -is cut off by the finite range of the interaction potential, whose precise effective value can be a function of the microscopic roughness. [33] In addition, the corresponding intensity decreases with decreasing h=a and thus disappears in the limit of h=a ! 0.

Asymptotic scaling
In many cases, nearly incompressible elastomers are treated as perfectly incompressible and their Poisson's ratio is approximated with ν ¼ 0:5. [4][5][6]35] However, analyzing the asymptotic behavior of � E for h=a ! 0 reveals a more concise picture, which is presented in the following.
The ratio � E=E � can only be a function of the two dimensionless numbers defining the problem, namely ν and h=a. The asymptotic limits for � E=E � at extreme and small confinement can be deduced from existing solutions for the considered confined elastomer. As will be shown in the remaining part of this Section 3.1, they turn out to be for h=a � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À 2ν p and ν < 0:5 3a 2 =ð32h 2 Þ for h=a � 1 and ν ¼ 0:5; where c 1 ðνÞ was introduced in Eq. (4). The condition ðh=aÞ 2 � 1 À 2ν for ν < 0:5 in Eq. (11) is motivated by the observation that for ν close to 0.5, the scaling of the function cðν; h=aÞ has two different small-h=a scaling regimes. The threshold between these regimes is characterized by the transition from a parabolic to a constant stress distribution as illustrated by Figure 1. Based on this we introduce the terminology that a film is extremely confined if h=a is small compared to ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À 2ν p and largely confined if this condition does not hold but h is still small compared to a.

Unconfined limit
The unconfined limit is nothing but a regular flat punch in contact with a semi-infinite half space. [14] Inserting its well-known stress-displacement relation F ¼ 2aE � u 0 into Eq. (7) yields the first case given in Eq. (11).

Extreme confinement limit
For h=a ! 0 and ν < 0:5, the proportionality between the elastic energy of a single height undulation and jũðqÞj 2 does not depend noticeably on the wave number q. Due to Parseval's theorem, a linear relation between local displacement and local stress occurs within most of the contact. This is why the stress is approximately constant within the contact in the bottom row of Figure 1 and zero outside so that its Fourier transform is Θð. . .Þ being the Heaviside step function, r ¼ ðx; yÞ the in-plane position with r ¼ jrj and J n ð. . .Þ the Bessel function of the first kind of order n. Using the stress-strain relation that follows from Eq. (1) yields a displacement at the origin u 0 ;uðr ¼ 0Þ of ) where h ¼ h=a is the reduced height and q ¼ q a. By definition, the l.h.s. of Eq. (15) and thus its r.h.s. is nothing but the inverse of the effective elastic modulus � E. The integral on the r.h.s. of Eq. (15) probably has no closed-form analytical solution and is even difficult to solve numerically because of the oscillations of the Bessel function. To reduce the effect of the oscillations, we rewrite the integral as f ðqÞ being the integrand on the r.h.s. of Eq. (15). In its rewritten form, the integral can be easily seen to have its dominant contribution from small q when h is small. This property does not change after using the appropriate small-q approximation for cðν;qÞ in the integrand. As a consequence, I can be solved analytically to be h =c 1 ðνÞ in the large-confinement limit for ν < 0:5, which ultimately translates into the corresponding expression stated in Eq. (11).

Large-confinement limit
For ν ¼ 0:5 and h sufficiently small compared to a, or, alternatively, ν < 0:5 and h=a in the intermediate scaling regime, the stiffness of a mode is no longer independent of q but instead proportional to q À 2 . As a consequence, the stress is no longer constant in the contact area but assumes the functional form [35] σðrÞ with σ 0 ¼ 2� σ, see also the second to last row in Figure 1. It should be noted that this stress distribution was originally derived for a finite elastomer of originally cylindrical shape sandwiched between two rigid planes and assuming a stick condition. However, Gent [35] already expected the functional dependence of � E on h=a for his set-up to be similar to that of films. Proceeding as in Section 3.1.2, we first determine σðqÞ to be The displacement in the origin then becomes As in the previous section, the integral cannot be solved analytically and the integrand oscillates too much to allow for a numerically robust integration. We therefore proceed again as described in the text around Eq. (16). This time, we did not identify a closed-form expression for the small-q expansion, but found a numerical value of 4 with six significant digits, so that we believe 4 to be the exact value for the integral on the r.h.s. of Eq. (20). Thus, comparing the r.h.s. of Eq. (20) with the definition of � E yields the large-confinement limit for � E and ν ¼ 0:5.
An exact representation of � Eðh=aÞ=E � may be achievable by adding (infinitely) many summands as they occur on the r.h.s. of Eq. (10). To make these sums satisfy the asymptotic limits, Eq. (22) must be generalized to sum rules. However, we did not find that proceedings along those lines appeared to be promising. Therefore, we will only use the single wave-number, asymptotically correct approximations for � Eðh=aÞ.

Deducing depinning force and range of stable crack growth from the effective modulus
Refs. [3][4][5][6]8] relate � Eðh=aÞ (or its inverse) to the energy release rate G, from which the pull-off force and crack propagation dynamics can be deduced. Similar to Yang and Li, [36] we start from the total energy formulation rather than the energy release rate, as we find this more direct and more intuitive.
The total potential energy of our system in an externally potential producing a constant force F reads where a negative value of F implies a (positive) tensile force causing a negative displacement u 0 . In this nomenclature, a c is the actual contact radius, which may now be different from the punch radius a.
In equilibrium, u 0 and contact radius a c both minimize the potential energy, i. e., @U=@a ¼ @U=@u 0 ¼ 0. In stable equilibrium, the Hessian produced by the second-order derivatives of U w.r.t. u 0 and a c must be positive definite.
When the system is displacement-driven, this condition reduces to @ 2 U=@a 2 c > 0. The generalized force acting on the radius a c , F a ; À @U=@a c , is easily deduced as This equation allows a c to be determined self-consistently for a given u 0 . However, it has to be kept in mind that a c cannot grow for positive F a when a c is equal to the punch radius a. This is why the case of a c ¼ a and a c < a must be treated separately. The contact radius starts shrinking when F a ðu 0 ; a c ¼ aÞ ¼ 0 À on retraction. Inserting this condition into Eq. (24) and solving for u 0 yields u 0 ¼ À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi which has to be evaluated at a c ¼ a to deduce the normal displacement at the point, where the contact is just about to start shrinking for the first time. The normal force acting on the punch can then be deduced from @U=@u 0 ¼ 0 to be Evaluating this force at a c ¼ a gives the depinning or pull-off force F p ; À FðaÞ, which is the maximum tensile force occurring right before the contact radius starts shrinking. Different representations of the same equation can also be found using the energy release rate and/or assuming the load-driven case. [2,16] As mentioned above, the previously determined contact radius a c is only stable if @ 2 U=@a 2 c > 0, which can be written in a convenient form that reproduces @U=@a c , which is 0 in equilibrium: The dividing line between stable and unstable crack propagation is defined by the condition @F a =@a c ¼ 0. It can be cast This criterion together with Eq. (21) can be used to explain why there is no stable crack growth in an adhesive contact between a flat punch and a semiinfinite elastomer. For confined bodies, especially when ν is close to 0.5, this procedure is no longer applicable, since the contact area is usually not circular. Figure 2 compares numerical results for � E=E � to the analytical results from the previous section using adjustable parameters α and β following from Eq. (22). It is shown that the simple scaling approach reproduces the overall trends fairly well. By design, the asymptotic limits are matched for ν < 0:5. Moreover, the location of the minimum in � Eðh=aÞ=E � , so it exists for a given Poisson's ratio, almost coincides between theory and simulation. However, the value of � E=E � in the minimum has an error of a few 10%. Errors are largest in the regime where a ν , < 0:5 elastomer shows similar behavior to an ideally incompressible solid. This can be rationalized by the h=a ! 0 asymptotics of a ν ¼ 0:5 body, which would require the parameters α and β to be redefined. For reasons of completeness, we note that the minimum of � Eðh=aÞ=E � for ν ¼ 0:5 is located at h=a ¼ 1:665 in our analytical treatment and at h=a � 1:23 � 0:02 in the GFMD data. The relatively large numerical uncertainty of the minimum location results from the minimum being shallow.

Effective modulus
To better resolve the discrepancies between the scaling approach and the numerical data, Figure 3 shows the ratio � E GFMD = � E s as a function of reduced height. For ν < 0:45, relative errors turn out to be quite insensitive to ν for any h=a. They can be approximated reasonably well with a single Gaussian constructed according to .
with A ¼ 0:46, B ¼ 0:75, and h 0 ¼ 2. A similar insensitivity of � E GFMD = � E s on ν holds for large ν only as long as h=a � 1. Interestingly, all � E GFMD = � E s curves almost coincide at � h=a ¼ 1, where they assume the value of 4=3 within a 3% margin. Unfortunately, the relative errors can exceed a factor of 1.5 for ðh=aÞ 2 > 1 À 2ν while h=a � 1. Nonetheless, they always remain below a factor of π 2 =2.
Since previous works [2,3,5,7,21] considered mostly stick conditions for the elastomer-punch interface, comparing our numerical data to existing data or (semi-) empirical approximations for � Eðh=aÞ may not appear meaningful at first sight. However, we note that the overall trends are similar and that comparisons between different boundary conditions and comparisons between FEM and GFMD data may yet be insightful. We find Hensel et al.'s [3] (well, ugly) fit function to match our data most beautifully, which is shown exemplarily for ν ¼ 0:495 in Figure 3, in particular for 0:5 < h=a < 10. Interestingly, for h=a < 0:5, their fit function is close to our analytical result. Thus, it does not capture the minimum associated with the large Poisson's ratios, which may well be because Hensel et al. assumed stick conditions between punch and elastomer.

Pull-off stress
In order to deduce the adhesive pull-off stress σ p ¼ F p =ðπa 2 Þ from � Eðh=aÞ, we replace @ ln � E=@ ln a c in Eq. (26) with À @ ln � E=@ lnh so that σ p ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 4 � EðhÞγ The derivatives are evaluated numerically from the GFMD data shown in Figure 2 using cubic spline interpolation. Results for the pull-off stress are shown in Figure 4. To ensure their correctness, the pull-off stresses were also computed for selected values of ν at h=a ¼ 0:1 with simulations mimicking tack tests. Agreement was always within 2%. Since direct adhesive simulations are much more demanding and more prone to discretization errors than computations of � E using non-overlap constraints, we believe the presented results to have errors well below 2%.
Using our analytical expression for � Eðh=aÞ directly to estimate the pull-off stress turned out to be disappointing. However, using the asymptotic scaling for σ p yields relatively satisfactory results, because it allows one to transition from ν < 0:5 to ν ¼ 0:5 scaling when crossing over from extreme to large confinement. Likewise, it is beneficial to transition from ν ¼ 0:5 to h=a ! 1 scaling when crossing over from confined to unconfined. The asymptotic solutions for σ p can be obtained from Eq. (30) by exploiting once more Eq. (11): σ p ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 8γE � =ðπaÞ p for h=a � 1 ð31aÞ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi c 1 ðνÞγE � = ah � � r for h=a � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À 2ν p and ν < 0:5 ð31bÞ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi for h=a � 1 and ν ¼ 0:5: ð31cÞ Related scaling relations have been proposed before [1,4,7,12,13,15,16,36] for varying boundary conditions (BCs), however, generally assuming ν ¼ 0:5. Historically first, Kendall [15] identified Eq. (31a) for the unconfined system using, as we do, a frictionless elastomer-punch interface. In the opposite limit, h=a ! 0, he found σ p / 1= ffi ffi ffi h p , which differs from our Eq. (31c) because Kendall used a slip condition for the elastomer-substrate interface, while we assumed a stick condition. Yang and Li [36] confirmed Kendall's scaling relation, albeit they corrected the numerical prefactor by multiplying Kendall's result with ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi , where K is the bulk modulus. In the case where Yang and Li employ our BCs, they also find Eq. (31c). In fact, Yang and Li considered all four possible combinations of elastomer-punch and elastomer-substrate BCs. However, they only considered ν ¼ 0:5.

Deducing ν from mechanical measurements
In order for our calculations to benefit the determination of the Poisson's ratio from mechanical measurements or, rather that of Δν ¼ 0:5 À ν, our results for � E and σ p are best represented as functions of ν for fixed values of h=a, as is done in Figure 5. This way E � and/or σ p only need to be determined once for a semi-infinite elastomer and once more for a confined elastomer. Similar approaches to determine Δν experimentally from � Eðh=aÞ have already been pursued successfully. [37,38] Information as that presented in Figure 5 is certainly only beneficial as long as we are not yet too deep in the extreme confinement limit since � E and σ p are no longer sensitive to logðΔνÞ in that regime.
Using very small h=a from the beginning is not necessarily effective either, since it might be equally important and infeasible to accurately align the flat punch as well as to account for the effects arising from the combined compliance of substrate, punch, and driving apparatus. Thus, for Δν suspected to exceed 10 À 3 and 10 À 4 , we would recommend to use h=a ¼ 0:05 and h=a ¼ 0:02, respectively.
Determining Poisson's ratios of confined layers to less than 10 À 4 might be possible through optical measurements of the bulge arising right next to the indenter. [39] However, we could not identify bulge characteristics, i.e., appropriately undimensionalized bulge widths or heights, which appear to be promising candidates. The thin slit between the indenter and elastomer seems to have the largest sensitivity to logðΔνÞ. Unfortunately, its determination would require extremely smooth surfaces and a high-accuracy measurement of the buried gap.

Crack formation and propagation
Crack growth during punch retraction becomes stable for sufficiently confined elastomers and large Poisson's ratios, i.e., below a critical film height h c ðνÞ. Evaluating the stability condition, Eq. (28), for ν ¼ 0:5, we locate the transition near h c ð0:5Þ=a ¼ 3:44 from the GFMD data and at h c ð0:5Þ=a ¼ 3:692 from our scaling ansatz. As ν decreases, the estimates for h c ðνÞ=a move to smaller values. However, h c ðνÞ=a becomes tedious to evaluate numerically from GFMD data once the elastomer is no longer very close to being incompressible.
Simulations [12,13] reveal elastic instabilities similar to those observed experimentally, thereby supporting the theoretical analysis. It yet seems unclear why and how the simulated patterns in two-dimensional contacts [13] break the symmetry of the mathematical problem, which has circular symmetry in the absence of discretization and periodic boundary conditions. To elucidate this issue further, we simulated the detachment process for three different Poisson's ratios using different ways in which the continuum model was terminated at small scales. Some of the most intriguing snapshots taken during detachment are compiled in Figure 6.
Every graph in a row in Figure 6 reflects the same continuum model in that ν and h=a is kept constant. However, they differ in terms of their discretization -leading to or suppressing lattice trapping at small scales -and in terms of the absence or presence of random noise, which is introduced with a Langevin thermostat. Despite representing the same continuum limit, all four graphs within a row look qualitatively different with the exception of the two right panels in the bottom row, which are both singly connected contact domains with the four-fold symmetry of the discretized model. The most highly symmetric patterns are obtained when thermal noise and lattice trapping are absent. While the ν ¼ 0:4 and ν ¼ 0:48 systems have circular symmetry, the ν ¼ 0:499 configuration reduces to a four-fold symmetry axis. The symmetry reduction is not due to the presence of periodic boundary conditions in our square domain but results from the discretization of the elastomer's surface into grid points forming a square lattice. We come to this conclusion because increasing the buffer between the punch and the boundary of the simulation cell does not change the point at which circular symmetry is broken. However, we observed that the rate of retraction can matter. For example, a complete loss of symmetry can occur even in the absence of thermal noise when decreasing from very small to extremely small retraction rates. Since our four-fold symmetry axis is only broken by the order in which numbers are added up, the complete loss of symmetry can only result from an accumulation of round-off errors. We suspect that a similar round-off error progression to significant digits is responsible for the low-symmetry configurations produced by Gonuguntla et al., [13] owing to them using a highly efficient conjugate gradient minimization method and/or because computers used smaller data precision in 2006 than they do nowadays. Given our results, we predict that instability patterns assume a quasi-circular symmetry, when elastomers are retracted quickly if the original surfaces are sufficiently planar.
Switching on temperature yield configurations similar to those observed experimentally and more so for a fine discretization avoiding lattice trapping. Specifically, the snapshots shown in the "no-trapping, T > 0" column resemble typical experimental images [4,[6][7][8][9] for ν > 0:45 and the contact shown for ν ¼ 0:4 in that column is similar to that depicted in Fig. 9a of Ref. [8] for ν ¼ 0:4. Due to lattice trapping, the non-contact patches show 90 � corners oriented w. r.t. the microscopic shape, similar to but substantially stronger than in the pioneering simulations by Gonuguntla et al., [13] who thus must have also discretized their domain into squares.
In the presence of thermal noise, the width of contact and of noncontact domains is similar in size. Their combined width indeed satisfies λ � 3h, which is the expected wavelength for nearly incompressible elastomers introduced at the beginning of this section. However, contact generally appears broader than non-contact due to (close-to) circular symmetry. This difference might matter for a comparison between singlewavelength pen-on-paper theory and real or realistic patterns. In addition, the surface tension can shift the characteristic wavelength to larger values. [9,13] We also analyze the effect of lattice trapping. For this phenomenon to occur, it does not matter whether the range of adhesion is decreased at fixed discretization Δx ¼ Δy or the mesh size is increased at fixed range of adhesion, as long as Δx is clearly less than typical contact and non-contact widths. Lattice trapping also counteracts symmetry reduction, as revealed most clearly in the bottom row of Figure 7, where thermal fluctuations are no longer strong enough to roughen the contact line during the course of the simulation. In other words, symmetry reduction can become an activated process in the case of lattice trapping whereby contact domains become thicker than non-contact regions upon retraction.
Unfortunately, our simulations do not correlate very well with some experiments regarding one aspect: for ν > 0:45, we usually observe the nucleation of non-contact below the punch center, while experiments often find fingershaped non-contact regions to emerge from the rim of the punch and then to move inward. [4,[6][7][8] Reasons for this discrepancy might be (i) the simulations ignore the effect of air pressure, which certainly favors the primary detachment to occur at the contact edge, (ii) the simulations neglect shear stress, which can be large near the contact periphery, (iii) no attempts were made to model viscoelastic effects, and (iv) the edge singularities in the normal stress are cut off too early due to a coarse discretization. Nonetheless, other experiments observed, as we did, crack nucleation in the center, [9,11,13,38] especially in cases where h=a � 0:1.
It is beyond the scope of this work to test all four hypotheses for why our non-contact domains nucleate in the centerin particular, as testing the first three does not fall into the realm of our model. However, we did investigate the fourth hypothesis by increasing the resolution from our default value to 4096 � 4096 for a h=a ¼ 0:1 punch and ν ¼ 0:495 elastomer while decreasing the range of adhesion and the rate of retraction so that the ratio of local elastic stiffness and maximum curvature of the tensile potential remained constant. As a consequence, the detachment nucleates at the periphery of the contact as is revealed in Figure 7 for the fine discretization, while it nucleated in the center for the coarser simulation.

Summary
In this work, we combine and streamline existing approaches to the mechanics of confined elastomers interacting with a flat punch. In doing so, we identify a relatively simple, yet physically motivated expression for how the effective modulus of the confined layer, � E, depends on its Poisson's ratio ν and the ratio of elastomer height and punch radius h=a. Using our approach, the Poisson's ratios no longer have to be close to 0.5. One consequence of this is that large confinement can be distinguished from extreme confinement, for which deviations from ideal incompressibility cannot be ignored. A central benefit of the pursued scaling ansatz is that the asymptotic dependence of � Eðh=a ! 0Þ allows the simulated pull-off force to be estimated reasonably well for any combination of ν and h=a. Of course, in real-laboratory experiments, the compliance of the substrate, the indenter, or, more generally speaking, the system must be considered when deducing � E in the extreme-confinement limit. Moreover, eliminating viscoelastic retardation implies (unrealistic?) requirements on the patience of experimentalists.
The central assumption of the analytical part of our study is that the elastic energy of a confined elastomer stems predominantly from surface undulations with wavelengths in the vicinity of the punch radius. This leads to a closedform expression in Eq. (10) for � Eðh=aÞ with two parameters of order unity, whose precise value can be fixed by demanding the asymptotic limits of h=a ! 0 and h=a ! 1 to be exactly reproduced. The pursued treatment can be repeated for boundary conditions (BCs) other than ours, which is a slip BC between the elastomer-punch interface and a stick BC for the elastomersubstrate interface. In these cases, the relation for the stiffness of the surface undulation of the wave vector q has to be derived or looked up in the literature, [11,[27][28][29]36] i.e., the replacement of Eq. (2). All remaining steps to estimate � Eðh=aÞ for other BCs can certainly be done by repeating the procedures worked out in this study. However, if both interfaces have slip boundary conditions, qualitatively different behavior ensues and other scaling relations apply than in the remaining three cases. [20,36] Investigating those is beyond the scope of this paper.
The analytical calculations are augmented with Green's function molecular dynamics (GFMD) simulations. They yield accurate reference data for the reduced elastic modulus and pull-off force as functions of reduced height and Poisson's ratios. In particular, the latter can be useful to determine experimentally the deviation of ν from ν ¼ 0:5.
The GFMD simulations also reveal that the initiation and the formation of cracks that occur during stable crack growth depend sensitively on the way in which continuum mechanics is terminated at small scales. For example, the (effective) range of interaction can determine whether cracks initiate from the punch center or from its periphery. Moreover, if interactions are so short ranged that lattice pinning ensues, small-scale features of the substrate, e.g. its crystallinity, can be reflected in the crack at coarse scales.