\section{Introduction}
 \subsection{Background}
The mechanics of a solid membrane is determined by a balance between its rigidity for in-plane (strain) and out-of-plane (bending) deformations.
For Graphene and other nanometer-thick crystalline two-dimensional (2D) membranes ({\emph{e.g.}} transition metal dichacogenides and black phosphorus), the in-plane stretching modulus $Y$ is very large, whereas the bending modulus $B$ is small, such that the characteristic length $ \ell_{bend} \approx \sqrt{\frac{B}{Y}}$ is much smaller than the system size $R_{\mathrm{sheet}}$ \cite{LWKH08,BBK11,QS13,CSZS14,RCCG15,Wetal16}.
Given the huge characteristic values of the von-Karman ratio, $vK =  (R_{\mathrm{sheet}}/\ell_{bend})^2$,
it is commonly assumed that
the bending rigidity does not affect the mechanics, and the response to exerted forces
is determined solely by the in-plane stiffness. While such an anticipation is justified when the exerted loads are purely tensile (\emph{e.g.} isotropic stretching of the sheet), it is obviously wrong to totally ignore the bending rigidity in the presence of compressive loads, as can be easily demonstrated by subjecting sheets to uniaxial compression  \cite{Milner89,Bowden98,Pocivavsek08,Huang10}. Here, the low bending rigidity underlies an instability of the compressed planar state, and the consequent formation of a strain-free buckled shape (if the sheet is suspended) or a wrinkle pattern (if the sheet is supported on a substrate) reflects the relevance of the
bending energy at scales much larger than $\ell_{bend}$. In this paper we study a conceptually similar, yet nontrivial
effect of the low bending rigidity in
indentation problems, where
radial tension induces compression in the azimuthal (hoop) direction, thereby making the weak bending energy a crucial player in the mechanical response of the sheet.
\\

Indentation experiments on suspended samples became a primary tool for measuring the stretching modulii of 2D materials \cite{NBK07,FTZM07,LWKH08,JP14,Fetal15,Letal15}. In a typical set-up, the sheet is supported on a thick, rigid substrate ({\emph{e.g.}} SiO), which contains a hole of radius $R \sim 1 \mu m$. A localized force is exerted by an AFM tip at the center of the suspended part of the sheet, and the force $F$ is measured as a function of the deflection $\delta$.
In most experiments \cite{LWKH08,Letal15} the stretching modulus, $Y$, is extracted by fitting the force-displacement curve, $F ( \delta )$, to a prediction of a ``membrane elasticity'' model, whereby the suspended
sheet is assumed to be clamped to the substrate at the edge of the hole.
This assumption implies that the indentation-induced stress field in the sheet is purely tensile,
and consequently has a dramatic influence on the estimated value of the stretching modulus \cite{Vella17}. However, an unequivocal, independent support for the validity of the clamping assumption has been lacking. Furthermore, since layers of Graphene (like graphite) are known to slide easily on each other (due to a very low inter-layer shear modulus), one may suspect that the interaction of graphene with a substrate is even weaker, such that the ``no-sliding'' assumption may not be satisfied.
 \subsection{Sliding, wrinkling, and response to applied forces.}
In order to understand the substantial effect of sliding on the indentation force, one must consider also the strength of the normal force that the substrate exerts on the supported part of the film. This interplay
can be demonstrated in a table-top example (Fig.~\ref{fig:demo}a):
attempting to push a tablecloth into a hole in a frictionless table, one finds that the tablecloth responds by changing its morphology - sliding towards the indenter, and forming radially-oriented blisters
(by buckling out from the table's plane) that release the hoop compression induced by the inward sliding.
Obviously, in such an experiment the fabric is not significantly stretched, indicating that
the combined effect of
in-plane sliding and out-of-plane deflection
may undermine
the use of indentation as a reliable probe for measuring the stretching modulus $Y$ of the sheet.
The theoretical model we introduce and analyze in this paper addresses the question that follows naturally from this simple observation:
{\emph{Under what conditions do sliding and deflection from the substrate
curb
the effect of the stretching modulus $Y$ on the indentation force
$F(\delta)$ ? }} \\

The mechanism for deflection from the substrate that we consider here, however, does not consist of blisters
(which are penalized by
surface energy and may be expected when the sheet-substrate attachment is sufficiently weak),
but rather of
small-amplitude wrinkles, such that the sheet-substrate distance $d$ remains close to its equilibrium value (see schematic Fig.~\ref{fig:schem-slide}d).
A central conclusion of our study is that when sliding and wrinkling are effective, the indentation force $F$ scales as:
\begin{align}
F &\sim \gamma_{\rm eff} \cdot \frac{R_{\mathrm{sheet}}}{R} \cdot \delta \ .
\label{eq:pseudo-linear-1}
\end{align}
Underlying Eq.~(\ref{eq:pseudo-linear-1}), which we call a {\emph{pseudo-linear}} response (and is valid above a certain threshold),
there is a highly non-linear geometric effect, comprising
a {\emph{global}} re-arrangement of the sheet in order to suppress the indentation-induced strain. The global nature of the pseudo-linear  response underlies its
dependence
on the overall size of the sheet, $(R_{\mathrm{sheet}})$,
in addition the hole's radius radius, $R$,
and an effective tension, $\gamma_{\rm eff}$, which is {\emph{independent}} on the stretching modulus $Y$, and may differ substantially from any pre-existing tension ($\sigma_0$) in the sheet. Specifically, $\gamma_{\rm eff}$ may reflect the bending rigidity and the steepness of the substrate-membrane VdW potential.
In contrast,
the standard linear response at infinitesimal indentation depth is
$F \sim (\sigma_0/R) \delta$ (up to logarithmic corrections \cite{Vella17}), being fully determined by the pre-tension $\sigma_0$, and the size of the suspended portion.

Before delving into the details of our model, let us provide a heuristic argument for the mechanism by which sliding and wrinkling give rise to a pseudo-linear response~(\ref{eq:pseudo-linear-1}). This argument is inspired by the example of indenting an ultrathin polymer sheet that is floating on a liquid bath \cite{Vella15,Paulsen16,Vella18b,Ripp20}.

\begin{figure}
\includegraphics[width=0.5\textwidth]{Fig1AD.png}
\includegraphics[width=0.5\textwidth]{Fig1B.png}
\includegraphics[width=0.5\textwidth]{Fig1C.png}
\caption{(a) Pushing a tablecloth into a hole does not cause significant stretching of the fabric, but rather its sliding towards the indenter, and the formation of a pattern of ``radial buckles''.
(b) Schematic sideview of indentation when a sheet is clamped to the hole's edge. The unavoidable stretching of radial lines
yields a tensile strain $\sim \sqrt{R^2 + \delta^2}/R - 1 \sim \delta^2/2R^2$. (c) If the sheet can slide inwards, the radial strain can be eliminated through
a displacement ${\rm u}_{\rm r} \!\sim\! -\delta^2/2R$ for $r>R$, and the ``bare" hoop compression thus acquired, $-{\rm u}_{\rm r}/r$, is relieved by undulations into the normal direction. (d) A schematic top view of a small portion of  a radially-wrinkled sheet. Since wrinkles cause a collapse of the hoop stress, force balance in the radial direction implies that $\sigma_{rr} \propto 1/r$; Boundary conditions at the far edge give: $\sigma_{rr} \approx \sigma_0 R_{\mathrm{sheet}} /r$.}
\label{fig:demo}
\end{figure}

\begin{figure}
\includegraphics[width=0.47\textwidth]{Fig2A.png}
\includegraphics[width=0.47\textwidth]{Fig2BC.png}
\includegraphics[width=0.47\textwidth]{Fig2D.png}
\caption{Schematic of our model system and key physical mechanisms. (a) A sideview of a solid membrane (sheet), which is supported on a smooth, rigid substrate in $R<r<R_{\mathrm{sheet}}$, and suspended in $r<R$. A pointwise indenter pushes at the center, causing the sheet to deflect downwards. Even though a sheet with finite bending modulus $B$ cannot accommodate a discontinuity of the tangent to its plane, our model allows the sheet to make a finite angle $\theta$ with the horizontal at the hole's edge. The reasoning is illustrated in panels (b) and (c), depicting two possible scenarios at the vicinity of the hole's edges (see App.~\ref{app:krr-BC}). In both scenarios, the tangent ``jumps'' over a short distance $\ell_{bc}^*$, Eq.~(\ref{eq:bendo-cap}), which constitutes a boundary layer of negligible energy cost. (d) We assume that the attractive force exerted by the substrate on the sheet is described by Zhang-Witten stiffness, Eq.~(\ref{eq:KsubWitten}), whereby the sheet-substrate distance remains in the VdW potential well. Delamination of the sheet from the substrate requires an energy barrier, and is not addressed within our model (see Sec.~\ref{sec:wrink-delam}).}
\label{fig:schem-slide}
\end{figure}

\subsection{Heuristic argument: stretching versus asymptotically isometric response \label{subsec:Heuristic}}
Let us contrast the two limit cases in the above example of pushing a tablecloth through a hole: {\emph{(a)}} perfect clamping of the sheet at the edge of the hole ($r=R$); or {\emph{(b)}} free sliding and ``wrinkling'' of the sheet on the substrate. \\

{\emph{(a) clamping:}} Assuming that prior to indentation the sheet is subjected to a uniform tension $\sigma_0$, the elastic energy associated with the work $F\cdot \delta$ of the indenter can be estimated as:
\begin{equation}
F \cdot \delta \sim U_{elas} \sim R^2 \left[\sigma_{0} + Y \left(\frac{\delta}{R}\right)^2 \right]  \left(\frac{\delta}{R}\right)^2  \ ,
\label{Eq:scaling1}
\end{equation}
where the stress in the sheet is estimated as the sum of the pre-tension $\sigma_{0}$ and the indentation-induced stress, $Y (\delta/R)^2$. Notice that the clamping assumption underlies our estimate of the radial strain, $\epsilon_{rr} \sim (\delta/R)^2$, as the indentation-induced extension of the radial distances (see Fig.~\ref{fig:demo}b). Note also that the bending energy is neglected, since we expect it to contribute only at some narrow, high-curvature zones, near the rim and around the indenter's tip. Equation~(\ref{Eq:scaling1}) shows that upon increasing $\delta$, the force transforms from a linear response, $F/\delta \sim \sigma_{0}$ (Column 2 of Table I), to a nonlinear, cubic response, $F / \delta^3  \sim Y/R^2$, which reflects the stretching modulus $Y$ \cite{Schwerin29} (Column 3 of Table I).
Actual calculations \cite{Vella17} yield a quantitative description of the transition between the two regimes (gray curve in Fig.~\ref{fig:response-1}).
\\

{\emph{(b) sliding and wrinkling:}}
Let us assume now that the sheet can slide freely on the substrate, such that material circles at radius $r$ undergo radial displacement $r \to r+ {\rm u}_{\rm r}(r)$. An inward displacement (${\rm u}_{\rm r}<0$) enables the sheet to retain the length of radial lines, thus avoiding the indentation-induced tensile strain $(\delta/R)^2$ in the radial direction; a simple calculation shows that
retention of the original length of radials, $R_{\mathrm{sheet}}$, requires a constant radial displacement
outside the hole ($r>R$):
\begin{equation}
{\rm u}_{\rm r}(r) \sim -\delta^2/R \
\label{eq:radial-scale}
\end{equation}
(Fig.~\ref{fig:demo}c).
Clearly, such an inward sliding causes a compression in the orthogonal planar (azimuthal) direction, since hoops of radius $r$ acquire a strain  ${\rm u}_{\rm r}/r$. If the normal attractive force exerted by the substrate is very strong, such a compression cannot be relieved, and the indentation force $F(\delta)$ is qualitatively similar to the clamping case discussed in the above paragraph. However, if the sheet can deflect even slightly from the substrate, then the compressive strain ${\rm u}_r/r$ can be eliminated by forming azimuthal undulations
whose characteristic wavelength may be very small, being determined by the bending modulus and the strength of sheet-substrate attachment. This scenario is the essential mechanism by which the table cloth in Fig.~\ref{fig:demo}a responds to the indentation force.

The elimination of tensile radial strain (by sliding) and compressive hoop strain (by deflection), suggests that the indentation force is not sensitive to the stretching modulus of the sheet.
Understanding this type of response, which involves only minute, asymptotically vanishing level of strain, and is thus called ``asymptotically isometric'' \cite{Vella15,Vella18,Paulsen2019,Davidovitch19}, is the essence of our manuscript. At a heuristic level, one can make progress by considering a small tension $\sigma_{0}$ pulling radially on the sheet at its far edge, $r=R_{\mathrm{sheet}}$, where $R_{\mathrm{sheet}} \gg R$. The presence of boundary tension implies that the stress in the sheet is not totally eliminated by sliding and wrinkling, and the response becomes dominated by the dependence of the residual stress on the indentation depth $\delta$. Since wrinkles eliminate the azimuthal component of the stress tensor, force balance on infinitesimal annular zones implies that there is a residual radial stress in the sheet $\sigma_{rr}(r) \approx \sigma_{0} R_{\mathrm{sheet}}/r$ (Fig.~\ref{fig:demo}d). The consequent divergence
at $r\to 0$ is resolved by the presence of an unwrinkled core of radius $L_I$ in which the stress saturates to its ``bare'' value $\sim Y\cdot (\delta/L_I)^2$. Continuity of radial stress at the boundary, $r=L_I$, between the wrinkled zone and the unwrinkled core, gives:
\begin{equation}
L_I \sim \frac{\sigma_{0}}{Y} \frac{R^2}{\delta^2}R_{\mathrm{sheet}} \ ,
\label{eq:Length-scale}
\end{equation}
(column 6 of Table II).

An interesting feature of such a sliding-wrinkling response is that only a negligible part of the indenter's work, $W_{inden} = F\cdot \delta$, is transmitted to the elastic energy of the sheet. In other words, the near-absence of strain, enabled by the combination of in-plane sliding
and out-of-plane deflection,
underlies a soft mode of an {\emph{asymptotically-isometric}} deformation, which eventually
controls the mechanical response to indentation.
This soft mode mechanics
may be realized by recalling that the only (finite) contribution to residual stress is the radial component $\sigma_{rr}$, whose integration yields $U_{elas} \sim \int_{L_I}^{R_{\mathrm{sheet}}} r \ dr \sigma_{rr}^2/Y \ \sim\  \sigma_0^2 R_{\mathrm{sheet}}^2/Y$, where we used Eq.~(\ref{eq:Length-scale}) and neglected logarithmic corrections and higher order terms associated with bending and other components of the stress. In contrast, the work done by the tensile load at the far edge against the indenter is $W_{ten} \sim 2\pi R_{\mathrm{sheet}} \sigma_{0} {\rm u}_r(R_{\mathrm{sheet}}) \sim \sigma_{0} \delta^2 R_{\mathrm{sheet}}/R $, where we used Eq~(\ref{eq:radial-scale}).
Introducing dimensionless parameters for the indentation depth and the sheet's radius:
\begin{equation}
\tilde{\delta} = \frac{\delta}{R} \sqrt{\frac{Y}{\sigma_0}} \ \  ; \ \ {\cal R} = \frac{R_{\mathrm{sheet}}}{R} \ ,
\label{eq:DG-00}
\end{equation}
we find that $W_{ten} > U_{elas}$ if:
\begin{equation}
\tilde{\delta} > \tilde{\delta}^{**}  \ \ \ {\rm where:} \ \ \tilde{\delta}^{**} \sim \sqrt{{\cal R}}
\label{eq:deltastst}
\end{equation}

 For $\tilde{\delta} >  \tilde{\delta}^{**}$, the elastic energy stored in the sheet can be neglected, and the indentation force is $F \approx \partial W_{ten}/\partial \delta$, so that we readily obtain Eq.~(\ref{eq:pseudo-linear-1}), with $\gamma_{\rm eff} = \sigma_0$.  \\

Crucially, Eq.~(\ref{eq:deltastst}) shows that
the depth $\delta$ required to reach such a pseudo-linear response
vanishes for a  ``nearly inextensible'' sheet ({\emph{i.e.}} $\sigma_0/Y \to 0$). In this regime, the combined effect of low energetic cost for bending and avoidance of indentation-induced strain, makes the solid sheet a ``bad capacitor'' of mechanical energy, and the work $W_{inden}=F\cdot \delta $ done by the indenter is transmitted almost entirely to the puller at the far edge \cite{Vella15,Vella18}:
\begin{eqnarray}
\tilde{\delta} < \tilde{\delta}^{**}:   \ \ \ F\cdot \delta  \ \ &\longrightarrow& \ \  W_{ten}  + U_{elas}    \nonumber \\
\tilde{\delta} > \tilde{\delta}^{**}:   \ \ \  F\cdot\delta  \ \ &\longrightarrow& \ \  W_{ten}  \ .
\end{eqnarray}

\subsection{Overview \label{subsec:overview}}
We introduce a minimal model to study the interplay between stretching, sliding, and wrinkling, and the
dependence of the indentation force
on actual physical parameters -- external tension, bending and stretching modulii of the sheet, and the strength of sheet-substrate attachment.
Asymptotic analysis of this model enables us to elucidate qualitatively different types of  response to indentation by identifying distinct parameter regimes.
\\

\subsubsection{Model and analysis}
Our model is depicted schematically in Fig.~\ref{fig:schem-slide}.
We consider a disk-like sheet of radius $R_{\mathrm{sheet}}$, with bending rigidity $B$ and Young modulus $Y$, which is suspended on a flat rigid substrate with a hole of radius $R \ll R_{\mathrm{sheet}}$ around it center ($r\!=\!0$), and a point-like indenter, which induces an out-of-plane deflection of amplitude $\delta$ at $r\!=\!0$.
We assume that the sheet is subjected to radial tension $\sigma_0$ at its far-edge $r = R_{\mathrm{sheet}}$, to which we will refer as ``pre-tension''. This may be an actual pre-tension ($\sigma_0=T_{\rm pre}$), applied prior to clamping the far-edge, or be exerted directly, such that the far edge, $r=R_{\mathrm{sheet}}$ is load-controlled rather than clamped. In our model, the normal force that resists deviations of the supported sheet from a planar state is characterized by a ``stiffness'' parameter $K_{\rm sub}$. Such a simplified response is known as Winkler foundation in the solid mechanics literature \cite{TimoshenkoBook}. The stiffness parameter $K_{\rm sub}$, together with the bending rigidity $B$ of the sheet, determine the deflections in the normal direction, which often take a periodic form that we call ``wrinkles'' -- the larger $K_{\rm sub}$ is, the smaller are the characteristic amplitude and wavelength of the emerging wrinkle pattern
\footnote{Familiar examples for $K$ are a liquid substrate (for a floating sheet), where deviations from planarity are resisted by the liquid gravity, hence $K_{\rm sub} = \rho_l g$  (where $\rho_l$ is the liquid density and $g$ is the gravity);
and a compliant substrate with a stiff near-surface layer of thickness $H$ and Young's modulus $E_s$, where $K \approx E_s/H$.}. For our primary interest here -- a highly rigid, undeformable substrate -- the stiffness $K_{\rm sub}$ was recognized by Zhang and Witten \cite{ZW07} as
\begin{equation}
{\rm Zhang\!-\!Witten \ stiffness:} \ \ K_{\rm sub} = V''(d_{\rm min}) \ ,
\label{eq:KsubWitten}
\end{equation}
where $V(d)$ is the attractive substrate-sheet potential per unit area, and $d_{\rm min}$ is the thermodynamic equilibrium distance between the sheet and the substrate (Fig.~\ref{fig:schem-slide}d). The Zhang-Witten stiffness assumes that the substrate is infinitely rigid, and the energetic cost for forming wrinkles (in addition to bending energy) is associated with the slight deviation of the sheet-substrate distance from its favorable value in the absence of any external loads. The assumption underlying this picture is that the VdW interaction is sufficiently strong, such that the energy barrier ({\emph{i.e.}} the depth of VdW potential well) that is necessary for the sheet to delaminate from the substrate cannot be reached. Instead, the small-amplitude undulations keep the sheet {\emph{everywhere}} within the VdW potential well of the substrate,
and the energetic penalty $\propto (d-d_{\rm min})^2$.
In Sec.~\ref{sec:wrink-delam} we will elaborate on the important difference between the relaxation of compression through such small-amplitude undulations and the formation of delamination zones, for which the energy cost per area is $\sim V(d_{min})$, independent of the actual sheet-substrate distance.

Our analysis is based on asymptotic analysis of the F\"oppl--von K\'arm\'an (FvK) equations, which describe the deformations of a thin solid sheet to exerted forces, assuming that the {\emph{local response}} is Hookean (namely, linear  stress-strain relationship), and that the deformed shape is characterized by small slopes. The FvK equations are {\emph{geometrically nonlinear}}, namely, the nonlinearity is {\emph{universal}} rather than material-dependent stemming from the coupling of out-of-plane deflections to in-plane strain.

In order to understand the response of the sheet to exerted loads, it is imperative to distinguish between the response to compressive and tensile stresses. If the stress exerted on a small piece of the sheet is (uniaxially or biaxially) tensile, the piece will stretch along the tension direction(s); we call this {\emph{tensile strain}}. In contrast, if the piece is under a compressive stress, it may buckle to reduce the compression level, and this mechanism gives rise to wrinkle patterns. If the sheet is sufficiently thin, or more precisely {\emph{highly bendable}}, the residual compression depends on the bending modulus and the exerted loads through a dimensionless parameter, called {\emph{bendability}} \cite{Davidovitch11}. The method by which the residual compression level is found, along with geometric features of the wrinkled state, has been called a \cite{Davidovitch11} {\emph{far from threshold}} analysis; this is an expansion of FvK equations around the singular limit of {\emph{tension field theory}}, which pertains to a compression-free sheet with no bending resistance ({\emph{i.e.}} $B=0$) \cite{Wagner35,MansfieldBook,Stein61,Pipkin86,Steigmann90}. It is thus crucial to understand that despite the smallness of the amplitude, the mere existence of wrinkles has a strong effect on the stress field in the sheet,
and therefore cannot be considered as a perturbation to some compressed, pre-buckled state.
\\

\begin{table*}
\centering
\caption{\label{Tab-I}
Summary of central results under various types of conditions of clamping and sliding at the hole's edge, upon increasing indentation depth $\delta$ (left to right), assuming the substrate stiffness $K_{\rm sub}$ is sufficiently large ($\beta \gg \tilde{\delta}^2$), such that the supported part of the sheet cannot wrinkle. The first row summarizes the response in the ``no-sliding'' case ($\sigma_{rr}(R) = \sigma_0$); the second row describes the ``free sliding'' case ($\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$), where compression develops above a threshold value, $\tilde{\delta}_c\approx 3.3$. The third row summarizes the effect of wrinkle formation in the suspended portion of the sheet.
}
\vspace{0.2cm}
\begin{ruledtabular}
\begin{tabular}{||c|c|c|c|c|c||}
&linear constant &cubic constant &asymptotic slope &asymptotic &comments \\
&$\frac{F}{\delta}$ & $\frac{F}{\delta^3}$ & at core's edge&tensile core   &  \\
& & & normalized by $\frac{\delta}{R}$ & & \\
\hline
no sliding &$- \frac{2 \pi}{\log ( \tilde{\delta} )} \sigma_0$ &$ \frac{0.166 \times 2 \pi}{R^2} Y$ &$ 0.63 $ &$L_I = R$ &pure tension  \\
& & & & &   (Sec.~\ref{subsec:II-clamping}) \\ \hline
sliding
&$- \frac{2 \pi}{\log ( \tilde{\delta} )} \sigma_0$ &$\frac{0.101 \times 2 \pi}{R^2} Y$ &$ 0.83 $ &$L_I \approx 0.6 R$ &wrinkling instability \\
no wrinkling& &  & & &$\tilde{\delta} \gtrsim 3.3$ (Sec.~\ref{subsec:II-sliding})  \\
\hline
sliding
& $- \frac{2 \pi}{\log ( \tilde{\delta} )} \sigma_0$  &$\frac{0.098 \times 2 \pi}{R^2} Y$ &$ 0.87$ &$L_I \approx 0.49 R$ & stable if $\beta \gg \tilde{\delta}^2$ \\
wrinkling  inside hole & & & & & (Sec.~\ref{subsec:II-wrinkling} )\\
 \end{tabular}
 \end{ruledtabular}
 \end{table*}
\begin{table*}
\caption{\label{Tab-II}
Summary of central results under various types of conditions of clamping and sliding at the hole's edge, upon increasing indentation depth $\delta$ (left to right), assuming $\beta <O(\tilde{\delta}^2)$, such that the substrate stiffness
$K_{\rm sub}$ is sufficiently small to allow relief of compression by forming wrinkles on the supported portion of the sheet.
The upper row describes the response under the ``sliding" BC ($\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$) in the parameter regime $\beta \ll 1$,
where the explicit values of the bending rigidity $B$ and substrate stiffness $K_{\rm sub}$ barely affect the stress field and  extent of the wrinkled zone. The middle row describes the response
when sliding is hindered by clamping the sheet at the far edge, such that ${\rm u}_{\rm r}(R_{\mathrm{sheet}}) = (1-\nu)\sigma_0 R_{\mathrm{sheet}}$. The bottom row describes the response under sliding (or hindered sliding) conditions, but when $1 \ll \beta \ll \tilde{\delta}^2$, such that the radial stress is governed by the residual hoop compression, $\sigma_{\theta\theta} \approx -2\sqrt{BK_{\rm sub}}$, rather than by the tensile load at the far edge. Here $\tilde{\delta} = \sqrt{Y/\gamma_{\rm eff}} \cdot R$.}
\vspace{0.2cm}
\begin{ruledtabular}
\begin{tabular}{||c|c|c|c|c|c|c||}
&linear &(sub) cubic & pseudo-linear & asymptotic slope &asymptotic &comments  \\
&constant & constant  & constant & at core's edge&tensile core   & \\
& $\frac{F}{\delta}$& $\frac{F}{\delta^3}$ & $\frac{F}{\delta}$ & normalized by $\frac{\delta}{R}$ & & \\
\hline \hline
sliding &$- \frac{ 2 \pi}{\log ( \tilde{\delta} )} \sigma_0$ & $- \frac{0.22 \times 2 \pi}{\log ( \tilde{\delta} )} \frac{Y}{R^2}$  & $\sim 2 \pi \sigma_0 {\cal R}$  & $1$
&$L_I \sim R \cdot {\cal R} \tilde{\delta}^{-2}$
& wrinkling instability  \\
$\beta \ll 1$ & & & & & & at $\tilde{\delta} \approx 3.3$ (Secs.~\ref{subsec:both-wrinkling},\ref{subsec:geo-limit})  \\ \hline
hindered sliding & $- \frac{2 \pi}{\log ( \tilde{\delta} )} \sigma_0$ & $\sim \frac{1}{\log {\cal R} } \frac{Y}{R^2}$ &none & $ 1 -O \left[ 1/\log{\cal R} \right] $ & $L_I \sim \frac{1}{\log {\cal R}} R$ & wrinkling instability  \\
$\beta \ll 1$
& & & & & & at $\tilde{\delta} \approx 3.3$ (Sec.~\ref{sec:Expansion}) \\  \hline
sliding & as above & as above & as above & as above  & as above &$\gamma_{\rm eff} = 2 \sqrt{B K_{\rm sub}}$ \\
$1 \ll \beta \ll \tilde{\delta}^2$
& $\sigma_0 \to \gamma_{\rm eff}$ & $\sigma_0\to \gamma_{\rm eff}$ & $\sigma_0\to \gamma_{\rm eff}$ & $\sigma_0\to \gamma_{\rm eff}$ & $\sigma_0\to \gamma_{\rm eff}$ &
(Sec.~\ref{sec:residual_compression})  \\
\hline
 \end{tabular}
 \end{ruledtabular}
 \end{table*}

\subsubsection{Classification of parameter regimes and central predictions}
Since the FvK equations are nonlinear, the stress cannot be considered a superposition of independent sources. Nevertheless, it is useful to identify three sources of stress that underlie the mechanical response.
\begin{equation}
\sigma_0 \ \ ; \ \ Y (\frac{\delta}{R})^2 \ \ ; \ \ 2\sqrt{K_{\rm sub} B}
\label{eq:stress-levels}
\end{equation}
The first source, $\sigma_0$, contributes a uniform isotropic tension to both radial and azimuthal (hoop) components of the stress tensor. The second source, $Y (\tfrac{\delta}{R})^2$, which is the only one that depends explicitly on the amplitude $\delta$,
gives rise to radial tension ({\emph{i.e.}} stretching radial lines) and hoop compression (pulling latitudes inwards). The last term,
$2\sqrt{K_{\rm sub} B}$, characterizes the residual hoop compression in the presence of radial wrinkles, namely, it is the minimal possible value, to which the hoop compressive stress can be suppressed with the aid of wrinkles
\cite{Cerda03,Davidovitch19}.
\\

The characteristic scales of stress (\ref{eq:stress-levels}) form
two dimensionless groups, in addition to $\tilde{\delta}$ and ${\cal R}$ (Eq.~\ref{eq:DG-00}), that we use to characterize the response to indentation at various parameter regimes:
\begin{equation}
\epsilon = \frac{B/R^2}{Y(\delta/R)^2} \ \ ; \ \ \beta = \frac{2\sqrt{BK_{\rm sub}}}{\sigma_0}  \ .
\label{eq:DG-0}
\end{equation}
The counterparts of the three dimensionless groups, $\tilde{\delta}, \epsilon, {\cal R}$,
have been used to describe the indentation of a floating ultrathin polymer sheet \cite{Vella15,Vella18}, whereas the additional parameter $\beta$, which describes the ratio between the residual compression and isotropic pre-tension, has received less attention in those studies.
The parameter $\tilde{\delta}$ -- a renormalized indentation depth -- is the ratio between the bare indentation-generated strain and the isotropic pre-tension in the sheet;
the parameter $\epsilon$ is the inverse of the ``geometric bendability'' -- the ratio between a (minimal) bending-related strain and the bare indentation-induced strain; the parameter ${\cal R}$ is the ratio between the lateral sizes of the sheet and the hole.

Throughout our study we will assume a highly bendable sheet, namely, $\epsilon \ll 1$, such that in-plane compression may be easily suppressed by wrinkling, and its size is large in comparison to the hole, namely ${\cal R} \gg 1$. Our primary interest is to understand the mechanics when $\tilde{\delta}$ is increased above a finite threshold value, $\tilde{\delta}_c \sim O(1)$, at which the indentation force is sufficiently strong to pull latitudes inwards and cause compression in part of the sheet. (Note that in the absence of pre-tension $\tilde{\delta}_c = 0$). In the rest of this introductory section we summarize the various types of response described by this model in terms of the parameters $\tilde{\delta}, \epsilon, \beta$, and ${\cal R}$. \\

\begin{center} {\emph{Regime (i)}} \end{center}
\begin{equation}
\beta \gg \tilde{\delta}^2
\label{eq:regime-i}
\end{equation}
In this parameter regime, the sheet-substrate attachment is so strong
that the supported portion of the sheet cannot relieve compression through wrinkling, even though the sheet is highly bendable.

If the sheet is not clamped to the hole's edge and can freely slide on the substrate, we find that for $\tilde{\delta}>\tilde{\delta}_c \approx 3.3$,  azimuthal (hoop) compression develops around $r=R$. For $\tilde{\delta}  \gg \tilde{\delta}_c$, the indentation-induced load dominates, and the stress becomes highly nonuniform and anisotropic, whereby the hoop-compressed zones extend upon increasing $\delta$. In the suspended part, the hoop compression can be effectively suppressed through the formation of radial wrinkles, but the supported part remains compressed.

The study of this parameter regime is the subject of the first parts of Sec.~\ref{sec:clamping_sliding} (A-C), and the results are summarized in Table I. Here, the central prediction of our study is a suppression of the force $F(\delta)$ due to sliding and wrinkling (second and third rows in Table I). Nonetheless, since the supported part of the sheet cannot wrinkle, the qualitative behavior -- transition from $F \sim \delta$ at $\tilde{\delta} \ll 1$ to $F\sim \delta^3$ at $\tilde{\delta} \gg 1$ -- is similar to the indentation of a sheet clamped at the hole's edge (first row in Table I). Note that the prefactor of the $\delta^3$ term changes significantly as function of the boundary conditions. Hence, even in this relatively simple regime, the extraction of the value of the Young modulus from an indentation experiment requires a careful consideration of the boundary conditions.
\\

\begin{center} {\emph{Regime (ii)}} \end{center}
\begin{equation}
\beta \ll 1
\label{eq:regime-ii}
\end{equation}
In this parameter regime, described
in Sec.~\ref{sec:clamping_sliding}.E-F,
the sheet-substrate attachment is sufficiently low such that it is energetically favorable to suppress hoop compression through radial wrinkles in both suspended and supported parts of the sheet. As a consequence,
the response to indentation
is qualitatively different from regime {\emph{(i)}}, and
is summarized in the first two rows of
Table II. A central prediction
is the emergence of a {\emph{pseudo-linear}} response, Eq.~(\ref{eq:pseudo-linear-1}), at sufficiently large indentation depth, $\tilde{\delta} \geq  O({\cal \sqrt{R}}) \gg 1$, which was motivated by our
heuristic discussion in Sec.~\ref{subsec:Heuristic}.
In Sec.~\ref{sec:Expansion} we discuss a situation where the sliding of the sheet on the substrate is hindered by clamping at the far edge, $r=R_{\mathrm{sheet}}$.
We find that far-edge clamping implies bi-axial tension
at the vicinity of $R_{\mathrm{sheet}}$,
even for large indentation depth ($\tilde{\delta} \gg {\cal R})$,
and thus eliminates the pseudo-linear response (second row of Table~II).
Nevertheless, the ability to relax hoop compression through wrinkling gives rise to dramatic suppression of the cubic response, $F\sim \delta^3$, when the clamping is at the sheet's edge ($r=R_{\mathrm{sheet}}$) in comparison to clamping
at the hole's edge ($r=R$). For $\tilde{\delta} \gg 1$, we find that the asymptotic ratio $F/\delta^3$ is proportional to $1/\log({\cal R})$.
\\

\begin{center} {\emph{Regime (iii)}} \end{center}
\begin{equation}
1 \ll \beta \ll \tilde{\delta}^2
\label{eq:regime-iii}
\end{equation}
In this parameter regime, the pre-tension $\sigma_0$ is irrelevant, and the substrate response is governed by a competition between the characteristic stress, $2\sqrt{BK_{\rm sub}}$, associated with the residual hoop compression in the wrinkled zone, and the bare indentation-induced stress, $Y(\tfrac{\delta}{R})^2$. We discuss this regime in Sec.~\ref{sec:residual_compression}, and show that the residual hoop compression gives rise to a comparable, {\emph{bending-induced}} radial tension \cite{hure12,Davidovitch19}.
This leads us to introduce an ``effective tension'' \cite{Tobasco20}:
\begin{equation}
\gamma_{\rm eff}
 \equiv \max\{\sigma_0,2\sqrt{BK_{\rm sub}}\} \ .
\label{eq:gamma-eff}
\end{equation}
Redefining the dimensionless amplitude:
\begin{equation}
 \tilde{\delta} = \frac{\delta}{R} \sqrt{\frac{Y}{\gamma_{\rm eff}}}  \ ,
 \label{eq:tdelta-2}
\end{equation}
we can characterize the mechanical response in regime {\emph{(iii)}} through a simple generalization of the predictions for regime {\emph{(ii)}, upon substituting (in all expressions that involve $\tilde{\delta}$): $\gamma_{\rm eff} = \sqrt{BK_{\rm sub}}$, rather than $\gamma_{\rm eff}=\sigma_0$. This is the content of the last row in Table II.
\\

\section{Clamping versus sliding and wrinkling \label{sec:clamping_sliding}}
We start by considering a perfectly axisymmetric deformation in response to indentation, namely, no wrinkles are allowed on the suspended or supported parts of the sheet. In Subsec.~\ref{subsec:II-clamping} and \ref{subsec:II-sliding} we address two types of boundary conditions (BCs). The first type is clamping at the hole's edge ($r=R$) with a ``pre-tension''
$\sigma_0$.
The second type of BCs allows for sliding of the sheet on the substrate, while a given tensile load, $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$, is exerted at the far edge of the sheet.
We show that the freedom to slide on the substrate significantly suppresses the indentation force.
Importantly, we find that if sliding is allowed, the sheet becomes azimuthally compressed in the vicinity of the hole's edge, if the dimensionless indentation depth, $\tilde{\delta}$, exceeds a critical value $\delta_c \approx 3.3$.
This indicates an instability to the formation of radial wrinkles, which we address in Subsec.~\ref{subsec:II-wrinkling}, assuming a sufficiently strong attachment to the substrate (regime {\emph{(i)}}, Eq.~\ref{eq:regime-i}), such that wrinkles can form only at the suspended part of the sheet. We use this case to introduce the basic principles of the {\emph{far from threshold}} approach \cite{Davidovitch11}, through which we characterize the emerging wrinkle pattern, and show how the formation of wrinkles underlies further, albeit modest suppression of the indentation force. In Subsec.~\ref{subsec:threshold-2} we relax the condition of infinitely strong sheet-substrate attachment, and find a second threshold, such that for $\tilde{\delta}>\tilde{\delta}^{*}(\beta)$, hoop compression is sufficiently strong to give rise to radial wrinkles also on the supported part of the sheet. In Subsecs.~\ref{subsec:both-wrinkling} and \ref{subsec:geo-limit} we address the parameter regime {\emph{(ii)}} (Eq.~\ref{eq:regime-ii}), where wrinkles expand throughout the supported part of the sheet and further suppress the indentation force, culminating with a transition to the pseudo-linear response, Eq.~(\ref{eq:pseudo-linear-1}).

\subsection{Clamping at the hole's edge}
\label{subsec:II-clamping}
Following Ref.~\cite{Vella17}, we address a circular sheet clamped with pre-tension $\sigma_0$ at the hole's edge, $r=R$.
The axial symmetry of the set-up calls for the use of polar coordinates. We denote the out-of-plane displacement by $z(r)$, and by $\psi(r)$ the radial derivative of the Airy stress function,
where the stress components are:
\begin{equation}
\sigma_{rr} = \frac{\psi}{r} \, , \, \sigma_{\theta\theta} = \psi'  \ .
\label{eq:define-Airy}
\end{equation}
The hoop component of the strain tensor $\varepsilon_{\theta\theta}$ and the consequent radial displacement field ${\rm u}_r$, satisfy:
\begin{gather}
\varepsilon_{\theta\theta} = \frac{{\rm u}_r}{r} = \frac{1}{Y} \left( \sigma_{\theta \theta} - \nu \sigma_{rr} \right) = \frac{1}{Y} \left( \psi' - \nu \frac{\psi}{r} \right) \nonumber \\
 \Rightarrow {\rm u}_r = \frac{1}{Y} \left( r \psi' - \nu \psi \right)
 \label{eq:disp-strain-stress-t}
 \end{gather}
whereas the radial strain is:
\begin{gather}
\varepsilon_{rr} = \frac{\partial {\rm u}_r}{\partial r}  + \frac{1}{2} (\frac{\partial z}{\partial r})^2 =
 \frac{1}{Y} \left( \sigma_{rr} - \nu \sigma_{\theta \theta}  \right) \nonumber \\
 = \frac{1}{Y} \left(\frac{\psi}{r} -\nu  \psi' \right) \ .
 \label{eq:disp-strain-stress-r}
\end{gather}
\\
The $2^{nd}$ FvK equation, expressing in-plane force balance (as well as compatibility of the stress and strain tensors with the displacement field), is:
\begin{align}
r \frac{d}{dr} \left[ \frac{1}{r} \frac{d}{dr} \left( r \psi \right) \right] &= - \frac{1}{2} Y \left( \frac{dz}{dr} \right)^2
\label{eq:FvK-dim-2}
\end{align}
and the $1^{st}$ FvK equation, which expresses force balance in the normal direction ($\approx \hat{z}$) is:
\begin{align}
\frac{1}{r} \frac{d}{dr} \left( \psi \frac{dz}{dr} \right) &= \frac{F}{2 \pi r} \delta ( r ) \ .
\label{eq:FvK-dim-1}
\end{align}
In the last equation we neglected a bending force, $B\partial^4z/\partial r^4$, due to the radial curvature of the sheet. As we explain in App.~\ref{app:krr-BC}, this term is significant only at the vicinity of the hole's edge, and its omission -- together with a suitable choice of BCs at the hole's edge -- is justified in all parameter regimes addressed in our paper (see schematic Fig.~\ref{fig:schem-slide}b-c).

Let us turn now to describe the BCs at the vicinity of the indenter ($r \to 0$), and the hole's edge:
\begin{align}
r \!\to\! 0  &:
(i) \,
z  \!=\!- \delta  \, \, \,
(ii) \,
{\rm u}_r \!=\! \lim_{r \rightarrow 0} \frac{1}{Y} \left( r \psi' \!-\! \nu \psi \right) \!=\! 0 \nonumber \\
r \!=\! R &:
(iii) \,
z \!=\! 0 \, \, \,
(iv) \, {\rm u}_r  \!=\! \frac{1}{Y}  \left(\! r \psi' \!-\! \nu \psi \!\right) \!=\! ( 1\! -\! \nu ) \frac{\sigma_0}{Y} R \ .
\label{eq:BC-clamping}
\end{align}
Note that clamping at $r = R$ in the presence of a pre-tension, $\sigma_0$, means that the radial displacement ${\rm u}_r (R )$  (rather than the load) is set to a fixed non-zero value, determined by the pre-indentation condition, as reflected in BC {\emph{(iv)}}, Eq.~(\ref{eq:BC-clamping}).
\\

Throughout this paper, we denote by $\psi , r , z$ dimensional values of the potential, radial (in-plane) length, and deflection (out-of-plane) length, respectively, and $\Psi , \rho$, and $\zeta$, for their dimensionless counterparts:
\begin{align}
\rho &= \frac{r}{R} , \, \, \, \Psi = \frac{\psi}{\sigma_0 R} , \, \, \, \zeta = \sqrt{\frac{Y}{\sigma_0}} \frac{z}{R} \ .
\label{eq:dimensionless-var}
\end{align}
Additionally, we define a dimensionless version of the force $F$:
\begin{align}
{\cal F} &
=  \frac{1}{2 \pi R} \sqrt{\frac{Y}{\sigma_0^3}}  \ F  \ .
 \label{eq:dimensionless-param}
\end{align}
The dimensionless form of the FvK equations (\ref{eq:FvK-dim-2},\ref{eq:FvK-dim-1}) is:
\begin{align}
\rho \frac{d}{d\rho} \left[ \frac{1}{\rho} \frac{d}{d\rho} \left( \rho \Psi \right) \right] &= - \frac{1}{2}  \left( \frac{d \zeta}{d \rho} \right)^2 \ ,
\label{eq:FvK-nondim-2}
\end{align}
\begin{align}
\frac{1}{\rho} \frac{d}{d \rho} \left( \Psi \frac{d \zeta}{d \rho} \right) &= \frac{\cal F}{2 \pi \rho} \delta ( \rho ) \ .
\label{eq:FvK-nondim-1}
\end{align}
The dimensionless version of the BCs (\ref{eq:BC-clamping}) is:
\begin{align}
\rho \rightarrow 0 \ &:  \, \, \, \,  (i ) \, \zeta = -\tilde{\delta} \, \, \, \, \, \, \, \, \,(ii ) \lim_{\rho\to 0} \left( \rho \Psi' - \nu \Psi \right) = 0 \nonumber \\
\rho =1 \ &:  \, \, \, \,  (iii) \, \zeta = 0 \, \, \, \, \, \, \, \, \, (iv) \, (\Psi' - \nu \Psi ) = 1 - \nu
\label{eq:BC-nondim-clamping}
\end{align}

Although Eqs.~(\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}) are nonlinear, there exists a transformation
\cite{Bhatia68}, which allows an analytic solution (up to integrals that can be evaluated numerically).
The analytic solution \cite{Vella17}, which we repeat in App.~\ref{app:axisymmetric}, enables us to express the force, $F$, the shape, $z(r)$, and the stress components, $\sigma_{rr}(r), \sigma_{\theta\theta}(r)$, for any value of $\tilde{\delta}$. These are shown, respectively, in the gray curves in Figs.~\ref{fig:response-1}-\ref{fig:shape}.  Let us discuss briefly some key features of of these results.
\\

At sufficiently large values of $\tilde{\delta}$, the force $F(\delta) \sim \delta^3$,
reflecting a transition from pretension-dominated stress (for $\tilde{\delta} \ll 1$)
to indentation-dominated stress: $\sim Y (\delta/R)^2$ (for $\tilde{\delta} \gg 1$).
The radial and hoop components of the stress, shown, respectively, by the dashed and solid curves in Fig.~\ref{fig:stress-profile-small} ($\delta = 3$) and Figs.~\ref{fig:hoop-stress-profile},\ref{fig:radial-stress-profile} ($\tilde{\delta} =10$),
indicate that an indented sheet clamped at the hole's edge is under pure tension at any indentation depth, in agreement with the discussion in
Subsec.~\ref{subsec:Heuristic}.

An interesting feature of the indentation force
is the absence of a true linear response (Fig.~\ref{fig:response-1}b). Instead, for $\tilde{\delta} \ll 1$ the response is sub-linear, with $F \sim -1/\log(\tilde{\delta})$. Such a sub-linear response appears also for the sliding BCs, albeit with a different numerical pre-factor (blue curve in Fig.~\ref{fig:response-1}). This peculiar feature emanates from the assumption of a point-wise indentation, and is intimately related to the  (integrable) divergence of the stress components at $r \to 0$ ($\sigma_{rr},\sigma_{\theta\theta} \sim r^{-1/3}$ \cite{Vella17}), which is observed in Fig.~\ref{fig:stress-profile-small}. For an indenter with a finite tip's radius $R_{tip}$, a linear response is recovered, with a numerical pre-factor that scales as $-[\log(R_{tip}/R)]^{-1}$ \cite{Vella17}.

Finally, it is noteworthy that the deformed shape (Fig.~\ref{fig:shape}) defers substantially from an ideal cone; this is signified by the slope in the vicinity of the hole's edge $(\tfrac{dz}{dr})_{r=R}$, which is only 63\% of the slope of an ideal cone (first row of Table I). We will show later that the slope at the vicinity of the edge depends strongly on the boundary conditions and other physical parameters, and may thus serve as an indirect experimental probe of the actual boundary conditions associated with a given set-up.

\begin{figure*}
\includegraphics[width=1.0\textwidth]{Fig3}
\caption{``Spring constant" $F/\delta$, {\emph{versus}} indentation amplitude $\delta$, where the force $F$ is normalized by $R^2\sigma_0$, and $\delta$ is normalized by $R\sqrt{\sigma_0/Y}$. Different colors represent different boundary conditions and physical parameters. Gray: clamping at the hole's edge. Blue: sliding of the sheet on the substrate, assuming the deformed shape is perfectly axisymmetric. Brown: sliding of the sheet on the substrate, where wrinkles are allowed to relax hoop compression only in the suspended part of the sheet (Regime {\emph{(i)}}, Eq.~\ref{eq:regime-i}). Green: sliding is hindered by clamping at the far edge, $r=R_{\mathrm{sheet}} \ (\approx 90 R$), and wrinkles are allowed to relax hoop compression in both suspended and supported parts of the sheet (Regime {\emph{(ii)}}, Eq.~\ref{eq:regime-ii}). Red: sliding throughout the whole sheet (same value of $R_{\mathrm{sheet}} \approx 90 R$), where wrinkles are allowed to relax compression in both suspended and supported parts of the sheet (Regime {\emph{(ii)}}, Eq.~\ref{eq:regime-ii}). For the problems that assume clamping (at the hole' edge (gray curve) or the sheet's edge (green curve)), $\sigma_0$ is the ``pretension" in the sheet, whereas for the sliding problems (blue, brown, and red curves), $\sigma_0$ is the tensile load exerted at the far edge. Regardless of the various BCs, the ``linear'' response at $\delta \ll R\sqrt{\sigma_0/Y}$, shown in panel B, is actually sub-linear, whereby $F/\delta \sim 1/|\log\delta| \to 0$ \cite{Vella17}. When wrinkles are not allowed on substrate (gray, blue, and brown curves) the asymptotic response at $\delta \gg R\sqrt{\sigma_0/Y}$ is cubic ($F/\delta\sim \delta^2$). When wrinkles can form on the substrate, the response become eventually ``pseudo-linear'' ($F/\delta \sim const$) after wrinkles can reach the sheet's edge (red curve), and sub-cubic ($F/\delta \sim \delta^2/|\log\delta|$) if the sheet's edge is clamped (green curve).
}\label{fig:response-1}
\end{figure*}
\begin{figure*}
\includegraphics[width=1.0\textwidth]{Fig4}
\caption{The radial stress at the hole's edge, $\sigma_{rr}(R)$, normalized by $\sigma_0$, {\emph{versus}} indentation amplitude $\delta$, normalized by $R\sqrt{\sigma_0/Y}$. Colors of curves represent different boundary conditions and physical parameters, as in Fig.~\ref{fig:response-1}.}
\label{fig:stress-edge}
\end{figure*}

\subsection{Sliding at the hole's edge and the buckling threshold in the suspended zone}
\label{subsec:II-sliding}
Now we address the axisymmetric (unwrinkled) state of the indented sheet in a set-up, where sliding of the sheet is allowed at the hole's edge (and on the substrate). Clearly, the only difference between this case and the above analysis of clamping at $r=R$, is encapsulated by the BC {\emph{(iv)}} in Eq.~(\ref{eq:BC-nondim-clamping}). For simplicity, we assume a fixed tensile load at the far edge, $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$. We note that, as long as $R_{\mathrm{sheet}} \gg R$, replacing this BC with clamping at the far edge (with pre-tension $\sigma_0$), gives rise to practically indistinguishable results.

In order to derive the appropriate BC at the hole's edge, we must consider the stress field in the supported part of the sheet. In this annular zone, $R<r<R_{\mathrm{sheet}}$, the sheet is subjected to radial tension $\sigma_0$ at $r=R_{\mathrm{sheet}}$ and an unknown radial tension $\sigma_{rr}(R)$ at $r=R$. This problem is readily recognized as the Lam\'e problem, and its classical solution yields the hoop and radial stress components \cite{TimoshenkoBook}:
\begin{equation}
\label{eq:Lame-stress}
\small{R\!<\!r\!<\!R_{\mathrm{sheet}}:} \left\{
\begin{array}{c}
  \sigma_{rr} (r) = \sigma_0 + \left(\sigma_{rr}(R) - \sigma_0\right)\frac{R^2}{r^2}     \\
   \sigma_{\theta\theta} (r) = \sigma_0 - \left(\sigma_{rr}(R) - \sigma_0\right)\frac{R^2}{r^2}
\end{array}
\right.
\end{equation}
(where we simplified the general solution for ${\cal R} = R_{\mathrm{sheet}}/R \gg 1$),
allowing one to express the radial displacement, ${\rm u}_r(r)$ as a function of $R_{\mathrm{sheet}},r,\sigma_0$, and $\sigma_{rr}(r)$:
\begin{equation}\
\!\!\!R<r<R_{\mathrm{sheet}}: \ {\rm u}_r(r) = \frac{r}{Y} [2 \sigma_0 - (1+\nu)\sigma_{rr}(r)]  \ .
\label{eq:dis-stress-Lame}
\end{equation}
Obviously, integrity of the sheet requires continuity of the radial displacement at the hole's edge, namely:
\begin{equation}{\rm u}_r(r \to R^+) = {\rm u}_r(r \to R^-) \ ,  \label{eq:matchingur}
\end{equation}
and the analogous relationship for the radial component of the stress reads:
\begin{equation}\sigma_{rr}(r \to R^+) = \sigma_{rr}(r \to R^-) \ . \label{eq:matchingsrr}
\end{equation}
In App.~\ref{app:krr-BC} we will elaborate further on the continuity of radial displacement and stress at the hole's edge and the validity of the corresponding equations (\ref{eq:matchingur},\ref{eq:matchingsrr}).

Equations~(\ref{eq:dis-stress-Lame},\ref{eq:matchingur},\ref{eq:matchingsrr}), together with Eqs.~(\ref{eq:define-Airy},\ref{eq:disp-strain-stress-t}) yield: $\psi/R + \psi' = 2\sigma_0$. Turning to dimensionless representation we obtain the BCs:
\begin{align}
\rho \rightarrow 0 \ &:  \, \, \, \,  (i ) \, \zeta = -\tilde{\delta} \, \, \, \, \, \, \, \, \,(ii ) \lim_{\rho\to 0} \left( \rho \Psi' - \nu \Psi \right) = 0 \nonumber \\
\rho =1 \ &:  \, \, \, \,  (iii) \, \zeta = 0 \, \, \, \, \, \, \, \, \, (iv) \, \Psi + \Psi' = 2
\label{eq:BC-nondim-sliding-axi}
\end{align}
The solution of the FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}), with the BCs~(\ref{eq:BC-nondim-sliding-axi}) can be obtained in a similar way to the solution in the preceding subsection (see App.~\ref{sub-app:sliding}),
allowing us to express the force, $F$, the shape, $z(r)$, and the stress components, $\sigma_{rr}(r), \sigma_{\theta\theta}(r)$, for any value of $\tilde{\delta}$. These are shown, respectively, in the blue curves in Figs.~\ref{fig:response-1}-\ref{fig:shape}.

One may notice that the qualitative behavior of the axisymmetric state with sliding BC's is very similar to the edge-clamped set-up. Considering the stress and force as functions of the dimensionless parameter $\tilde{\delta}$ (Figs.~\ref{fig:response-1},\ref{fig:stress-edge}, respectively), their magnitudes scale similarly in both set-ups, for $\tilde{\delta} \ll 1$ as well as for $\tilde{\delta} \gg 1$, whereas the numerical pre-factors become smaller once sliding is allowed. Intuitively, sliding allows the sheet to moderately relax the stretching in the suspended part at the expense of more stretching at the supported part. An interesting observation is the pronounced
effect of the change in BCs on the shape (Fig.~\ref{fig:shape}). Specifically, we found that for $\tilde{\delta} \gg 1$, the contact angle of the membrane at the edge approaches the asymptotic value: $(\tfrac{dz}{dr})_{r=R} \rightarrow 0.82 \frac{\delta}{R}$ -- an increase of over 25\% of the slope under clamped BC's.

A dramatic feature of the stress profile is the emergence of an annular zone on both sides of the edge, where the hoop stress is compressive (blue curve in Fig.~\ref{fig:hoop-stress-profile}), for $\tilde{\delta} \gtrsim 3.3$.
While a compressive stress may not be relieved in the supported part due to a strong attachment to the substrate (regime {\emph{(i)}}, Eq.~\ref{eq:regime-i}), the existence of compression in the suspended part of a thin sheet clearly gives rise to a wrinkling instability. Understanding the wrinkle pattern under such physical conditions is the subject of the next subsection.

\begin{figure}
\includegraphics[width=0.5\textwidth]{Fig5.pdf}
\caption{(a) The hoop stress (solid) and radial stress (dashed) for a sheet that is clamped (gray) or slide (blue) at the hole's edge. Here, the value of the dimensionless indentation depth $\tilde{\delta} = 3$, for which both cases the sheet is under pure tension, and the axisymmetric response is stable. Both stress components are normalized by a constant $\sigma_0$ (see text). (b) The profile of the suspended sheet. Here, radial distances are normalized by the hole's radius $R$, whereas vertical distances are normalized by the indentation depth $\delta$. Note the deviations in both cases from a perfect conical shape.}
\label{fig:stress-profile-small}
\end{figure}

\subsection{Wrinkling in the suspended zone}
\label{subsec:II-wrinkling}
We come to study the simplest case in which wrinkles affect the mechanical response, where the hoop compression induced by sliding at the hole's edge gives rise to wrinkles in the suspended portion of the sheet ($r<R$), but not at the supported part ($r>R$).
In order to study the effect of wrinkles on the stress and thereby on the indentation force, we employ {\emph{tension field theory (TFT)}} \cite{Wagner35,MansfieldBook,Stein61,Pipkin86,Steigmann90}. In this approach, one assumes that wrinkles suppress almost entirely compressive stress, such that one of the two principal components of the
stress tensor in the wrinkled zone is positive,
corresponding to tensile stress {\emph{along}} wrinkles, whereas the other principal component vanishes, signifying the direction along which wrinkles undulate.
The stress field in the whole sheet is then obtained by matching the displacement field and the compression-free stress in the wrinkled zone to the adjacent, purely tensile zones, where both principal stress components are non-negative.

Applying the TFT methodology to our indentation problem, the axial symmetry of the set-up suggests that for $\tilde{\delta} > 3.3$, confinement of latitudes occurs in an annular zone, $L_I<r<L_{\mathrm{out}}$, where $L_I <R$, and $L_{\mathrm{out}} >R$. In this subsection we assume that in the supported part, $R<r<L_{\mathrm{out}}$, the large effective stiffness $K_{\rm sub}$, Eq.~(\ref{eq:KsubWitten}),
prohibits the formation of wrinkles, such that the sheet must accommodate the indentation-induced hoop compression; however, in the suspended part, $L_I<r<R$, the formation of radial wrinkles underlies {\emph{collapse}} of hoop compression. Hence, the sheet is naturally divided into three parts: {\emph{(i)}} $R<r<R_{\mathrm{sheet}}$ -- where the supported sheet undergoes a planar axisymmetric deformation;
{\emph{(ii)}} $L_I<r<R$ -- where the suspended sheet is wrinkled; {\emph{(iii)}} $r<L_I$ -- where the suspended sheet is unwrinkled and the stress is purely tensile. In the sequel, we will describe the stress and deformation in each zone
and the matching among them.
\\

\noindent
\underline{Zone {\emph{(i)}} $R<r<R_{\mathrm{sheet}}$:}
Similarly to Subsec.~\ref{subsec:II-sliding}, the state of the sheet in this part is determined by solving the planar Lam\'e problem, subject to radial tensile load at the far edge, $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$, and a radial tension $\sigma_{rr}(R)$, which must be determined by matching the three zones. In this zone, the stress components are given by Eq.~(\ref{eq:Lame-stress}) through the unknown $\sigma_{rr}(R)$, and the ratio between the radial displacement and stress is given by Eq.~(\ref{eq:dis-stress-Lame}), which we repeat here for completeness:
\begin{equation}\
{\rm u}_r(r) = \frac{r}{Y} [2 \sigma_0 - (1+\nu)\sigma_{rr}(r)]    \ .
\label{eq:dis-stress-Lame-1}
\end{equation}

\noindent
\underline{Zone {\emph{(ii)}} $L_I<r<R$:} Here, the formation of wrinkles underlies a collapse of the hoop compression, such that we need to solve the radial force balance equation with $\sigma_{\theta\theta}=0$. (More precisely, a TFT solution is the leading order in a ``high bendability'' expansion, $\epsilon \to 0$, of the FvK equations \cite{Davidovitch11},
rather than a standard expansion around the compressed, axisymmetric state \cite{TimoshenkoBook}). Technically,  in the wrinkled zone radial force balance is obtained by satisfying Eq.~(\ref{eq:define-Airy}) with $\psi(r) = constant$, whereas the relationship $\varepsilon_{\theta\theta} = {\rm u}_{\rm r}/r$ (\ref{eq:disp-strain-stress-t}) between the hoop strain and radial displacement is ``ignored''  since it merely determines a comparable contribution to the hoop strain ($\varepsilon_{\theta\theta} = -\nu \varepsilon_{rr}$), which is missing from the RHS of Eq.~(\ref{eq:disp-strain-stress-t}) due to the excess length in the wrinkly undulations \cite{Davidovitch11}. Equation (\ref{eq:FvK-dim-2}) also relies on (\ref{eq:disp-strain-stress-t}) and its validity is thus limited to an axisymmetric, unwrinkled state, hence it is likewise ignored. Requiring continuity of radial displacement and stress, Eqs.~(\ref{eq:matchingur},\ref{eq:matchingsrr}),
as well as continuity of the deflection $z(R) = 0$, and employing Eqs.~(\ref{eq:define-Airy},\ref{eq:disp-strain-stress-r},\ref{eq:FvK-dim-1}), we obtain the stress components, deflection $z(r)$ and radial displacement in the wrinkled zone:
 \begin{eqnarray}
 \sigma_{rr}(r) &=& \sigma_{rr}(R)\frac{R}{r} \, , \,
 \sigma_{\theta\theta}(r) = 0 \ \Rightarrow \ \psi(r) = R \cdot \sigma_{rr}(R)
 \nonumber \\
 z(r) &=& \arctan\theta \cdot ( r - R ) \approx \theta \cdot ( r - R )  \ ,  \nonumber  \\
 {\rm u}_{\rm r}(r) &=& \!- \!\frac{1}{2} \theta^2\cdot (r\!-\!L_I) \!+\! \frac{1}{Y}{\sigma_{rr}(R)} R \log \left(\! \frac{r}{L_I} \!\right)
\nonumber \\
&+&  {\rm u}_{\rm r}(L_I)
 \label{eq:wrinkled-I}
 \end{eqnarray}
where $\theta \ll 1$ is the angle between the suspended sheet and the planar substrate at $r=R$, and ${\rm u}_{\rm r}(L_I)$ is the radial displacement at the
edge of the wrinkled zone, $r=L_I$. We re-emphasize that although both radial displacement and Airy potential in the wrinkled zone are given by axisymmetric functions, the presence of symmetry-breaking wrinkles is reflected in the violation of the relationship (\ref{eq:disp-strain-stress-t}) between , ${\rm u}_{\rm r}(r)$ and $\psi(r)$.
\vspace{0.2cm}

\noindent
\underline{Zone {\emph{(iii)}} $0<r<L_I$:}
In the purely tensile core the state is again axisymmetric, and the FvK equations there, expressed through the dimensionless functions
$\zeta(\rho)$ and $\Psi(\rho)$~(Eq.~\ref{eq:dimensionless-var}) are correspondingly given by Eqs.~(\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}), with the strain-displacement relationship for both parts of the strain tensor (\ref{eq:disp-strain-stress-t},\ref{eq:disp-strain-stress-r}).
Exploiting once again the continuity of the radial stress $\sigma_{rr}(r) = \psi(r)/r$ and the deflection $z(r)$,
we obtain the BCs:
\begin{eqnarray}
\zeta(0)  = - \tilde{\delta}  \ \ \ &;& \  \ \
\lim_{\rho \to 0}\left( \rho \Psi' - \nu \Psi\right) = 0 \nonumber \\
\zeta(\tilde{L}_{I}) =\tilde{a} (\tilde{L}_{I}-1) \ \ \ &;& \ \   \Psi(\tilde{L}_{I}) = \Psi(1) = \frac{\sigma_{rr}(R)}{\sigma_0}
\label{eq:BC-sliding-wrinkling-I}
\end{eqnarray}
where $\tilde{L}_{I}$ and $\tilde{a}$ are dimensionless versions of the core radius and the slope at the hole's edge:
\begin{equation}
\tilde{L}_{I} = L_I/R \ \ ; \ \ \tilde{a} = \sqrt{Y/\sigma_0} \theta \ .
\label{eq:new-nondim-param}
\end{equation}
Similarly to the previous subsections, we find that the nonlinear Eqs.~(\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}) with the BCs~(\ref{eq:BC-sliding-wrinkling-I}) can be solved analytically. Namely, for a given value of the control parameter $\tilde{\delta}$ and given values of the three unknowns, $\Psi(1), \theta$, and $\tilde{L}_{I}$, there is a single analytic solution that fully characterizes the function $\Psi(\rho)$, the related components of the stress (\ref{eq:define-Airy}), and the deflection $\zeta(\rho)$, in the interval $0<\rho<\tilde{L}_{I}$. Since the state in this core zone is axisymmetric (unwrinkled), Eq.~(\ref{eq:disp-strain-stress-t}) implies that the radial displacement at $r=L$ satisfies:
\begin{gather}
{\rm u}_r(L_I) = \frac{1}{Y} \left( L \psi'(L_I) - \nu \psi(L_I) \right)  \nonumber \\
= \frac{\sigma_0}{Y} R \left( \tilde{L}_{I} \Psi'(\tilde{L}_{I}) - \nu \Psi(\tilde{L}_{I}) \right) \ .
 \label{eq:ur-core}
\end{gather}

\noindent
\underline{Matching conditions:}
In addressing the zones {\emph{(i-iii)}}, we only used the continuity of radial stress (hence $\Psi(\rho)$), and the deflection $\zeta(\rho)$. In order to determine the three unknown variables, $ \Psi(1)
\theta$, and $\tilde{L}_{I}$, we must invoke three other matching conditions. Two of them are continuity of the slope, $\zeta'(\rho)$, and hoop stress, $\sigma_{\theta\theta}(r) = \psi'(r) = \sigma_0 \Psi'(\rho)$, at the borderline between the tensile core and the wrinkled zone,
yielding two equations:
\begin{equation}
\zeta'(\tilde{L}_{I}) = \tilde{a} \ \  \ \ ; \  \ \ \ \Psi'(\tilde{L}_{I}) = 0 \ .
\label{eq:twomore}
\end{equation}
(As was noted in the similar problem of indenting a floating sheet \cite{Vella15}, these two equations do not follow from local force balance at $r=L_I$ {\emph{per se}}, but rather from minimization of the total energy of a wrinkled state, which is realized when the hoop stress is continuous throughout the sheet). The last matching condition is the continuity of radial displacement at the hole's edge, $\rho=1$, which is obtained through Eqs.~(\ref{eq:dis-stress-Lame-1},\ref{eq:wrinkled-I},\ref{eq:ur-core}, \ref{eq:twomore}), yielding:
\begin{equation}
\Psi(1) \cdot (1-\log\tilde{L}_{I})  - \frac{1}{2}\tilde{a}^2\cdot (1-\tilde{L}_{I})  = 2  \ .
\label{eq:thirdmatching}
\end{equation}
\\
With the three equations~(\ref{eq:twomore},\ref{eq:thirdmatching}), and the four BCs~(\ref{eq:BC-sliding-wrinkling-I}), the FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}), which are two coupled $2^{nd}$ order ODEs, yield a single solution for
$\Psi(\rho), \zeta(\rho)$ in the interval $0<\rho<\tilde{L}_{I}$, as well as the three unknowns, $\tilde{L}_{I}, \tilde{a}, \Psi(1)$. The details of the analytic solution are given in App.~\ref{sub-app:wrinkling}. Together with Eqs.~(\ref{eq:Lame-stress},\ref{eq:dis-stress-Lame-1},\ref{eq:wrinkled-I},\ref{eq:new-nondim-param}), this solution fully characterizes the displacement and stress fields
for any value of the dimensionless control parameter $\tilde{\delta} >\tilde{\delta}_c$.

The brown curves in Figs.~\ref{fig:response-1}-\ref{fig:stress-edge} show the force, $F$, and the radial stress at the hole's edge, $\sigma_{rr}(R)$, upon increasing $\delta$ above $\tilde{\delta}_c \approx 3.3$, and the brown curved in Figs.~\ref{fig:hoop-stress-profile}-\ref{fig:shape} show the hoop and radial stresses, $\sigma_{\theta\theta}(r),\sigma_{rr}(r)$, and the
shape, $z(r)$, at the suspended part,
for $\tilde{\delta}=10$.
One may notice that for any $\tilde{\delta} > 3.3$, wrinkling in the suspended part of the sheet reduces slightly further the force (in comparison to the unstable, axisymmetric deformation with sliding, represented by the blue curve).
Note also that the formation of wrinkles acts to slightly increase the angle $\theta$ at the hole's edge
in comparison to an unwrinkled deformation and
to extend the azimuthally-confined zone (fourth and fifth columns in Table 1, respectively).
Intuitively, since wrinkles suppress the energetic cost of hoop strain, it is favorable to extend this zone.

One should note the discontinuity exhibited in Fig.~\ref{fig:hoop-stress-profile} by the hoop stress at the hole's edge. Such a discontinuity does not violate any force balance, and is therefore physically allowed. More precisely, while this discontinuity emerges naturally in TFT, which describes the infinite bendability limit (of a hypothetical sheet with no bending rigidity, {\emph{i.e.}} $\epsilon \!=\! 0$),  we do expect the formation of a  ``boundary layer'' at the vicinity of the hole's edge, whose length vanishes as $\epsilon \to 0$, over which the ``jump'' in the hoop stress occurs (similarly, but not identically, to the boundary layer that regularizes a ``jump'' in the radial stress, see App.~\ref{app:krr-BC}). Nevertheless, the consequent effect on the elastic energy is negligible, {\emph{i.e.}} ``sub-dominant'' in the terminology of the far-from-threshold approach \cite{Davidovitch11}.

Taken together, these results demonstrate the wrinkles-assisted
suppression of the energetic cost of the deformation induced by indentation, and consequently a reduction of the force constant $F ( \delta ) / \delta^3$ in comparison to
the analogous nonlinear force constant for the unwrinkled state. These observations reinforce our qualitative discussion in Sec.\ref{subsec:Heuristic}, indicating that the formation of wrinkles implies a non-perturbative modification to the stress field, and thereby to the indentation force.

\begin{figure*}
\includegraphics[width=1.0\textwidth]{Fig6.pdf}
\caption{The hoop stress
for a dimensionless indentation depth $\tilde{\delta} = 10$ (left-suspended part, right - supported part). Distances are normalized by the hole's radius $R$ and stress is normalized by $\sigma_0$ (see text). The colors correspond to the various types of BCs, noted already in the caption of Fig.~\ref{fig:response-1}: gray (clamping at the hole's edge); blue (axisymmetric (unstable) response under sliding at the hole's edge); brown (wrinkling at the suspended part of the sheet only); green (wrinkling in both suspended and supported parts of the sheet, for a sheet with ${\cal R} \approx 90$ ($\tilde{\delta}^{**} ({\cal R}=90)<10$, Eq.~(\ref{eq:deltastst}));
red (wrinkling in both suspended and supported parts of the sheet, for a sheet with ${\cal R} = 6$ ($\tilde{\delta}^{**} ({\cal R}=6)>10$)). Note that only the brown curve is discontinuous at the hole's edge.}
\label{fig:hoop-stress-profile}
\end{figure*}

\begin{figure*}
\includegraphics[width=1.0\textwidth]{Fig7.pdf}
\caption{Same as Fig.~\ref{fig:hoop-stress-profile}, but for the radial stress. Note that all curves are continuous at the hole's edge.}
\label{fig:radial-stress-profile}
\end{figure*}

\subsection{Buckling threshold in the supported zone \label{subsec:threshold-2}}
In the previous subsection we let wrinkles suppress hoop compression only in the suspended part of the sheet,
whereas the supported part of the sheet remains unwrinkled.
In order to identify the parameter regime at which such a scenario may be realized, we note that the supported sheet is subjected to hoop compression at the vicinity of the hole's edge that keeps increasing in magnitude and spatial extent upon increasing the indentation depth. Physically, such a state is mechanically stable if the hoop compression is below the threshold value, $\approx 2\sqrt{BK_{\rm sub}} = \beta \sigma_0$ (Eq.~\ref{eq:DG-0}),
at which the supported sheet buckles.
This criterion is well known for uniaxial deformations \cite{Milner89,Bowden98,Pocivavsek08,Huang10} and was shown to be relevant also for more complicated, non-uniaxial confinement problems \cite{Davidovitch19,Bella17}). Considering our solution in Subsec.~\ref{subsec:II-wrinkling}, we note that the hoop compression at the edge ($\sigma_{\theta\theta}(R) = \psi'(r \to R^+) $) is approximately $ 0.11 \cdot Y(\delta/R)^2 = 0.11 \sigma_0 \tilde{\delta}^2$ (where we assumed $\tilde{\delta} \gg 1$ for simplicity). Hence, we obtain that the indentation depth, $\tilde{\delta}^{*}(\beta)$, at which the {\emph{supported}} part of the sheet becomes wrinkled is given by:
\begin{equation}
\tilde{\delta}^{*} (\beta) \approx \left\{ \begin{array}{cc}
 \sqrt{0.11}\cdot  \beta^{1/2}    &   \beta \gg 1 \\
    \tilde{\delta}_c \approx 3.3 &  \beta \ll 1
\end{array} \right.
\label{eq:threshold-stst}
\end{equation}
Note that for $\beta <1$, the resistance to buckling in the supported part is sufficiently low, such that both supported and suspended parts of the sheets become wrinkled almost simultaneously, as soon as indentation-induced hoop compression emerges at $\tilde{\delta}
\gtrsim \tilde{\delta}_c \approx 3.3$.
\\

Equation~(\ref{eq:threshold-stst}) shows that the analysis in Subsec.~\ref{subsec:II-wrinkling} describes the parameter regime $\beta \gg 1 \ \& \  \tilde{\delta}_c < \tilde{\delta} \ll \tilde{\delta}^{*} (\beta)$, namely, where the sheet-substrate attachment is sufficiently strong to prevent wrinkling in the supported part, for sufficiently small indentation depth. This is precisely regime {\emph{(i)}} we described in Subsec.~\ref{subsec:overview}.

In the rest of this section, we will turn our attention to
regime {\emph{(ii)}}, $\beta \ll 1 \ \& \  \tilde{\delta} > \tilde{\delta}_c$, at which both suspended and supported parts of the sheet become wrinkled at $\tilde{\delta} \gtrsim \tilde{\delta}_c$, and the sheet-substrate attachment
does not affect the residual stress field.
In Sec.~\ref{sec:residual_compression} we will address regime {\emph{(iii)}}, $\beta \gg 1 \ \& \  \tilde{\delta} \gg \tilde{\delta}^{*}(\beta)$, at which the residual compression in the wrinkled, supported part of the indented sheet must be taken into consideration.

\subsection{Wrinkling in both suspended and supported zones \label{subsec:both-wrinkling}}
Considering the parameter regime {\emph{(ii)}}, $\beta\ll 1$ and $\tilde{\delta} > \tilde{\delta}_c \approx 3.3$,
we follow our analysis in Subsec.~{\ref{subsec:II-wrinkling}}, noting
that since $\beta \ll 1$, the direct effect of the sheet-substrate attachment on the stress field in the sheet is negligible, and therefore
the standard TFT approach of Subsec.~{\ref{subsec:II-wrinkling}} can be employed also here. Namely -- in the wrinkled zone, $L_I<r<L_O$, the stress field is given by a tensile radial stress, $\sigma_{rr}(r)>0$, and negligible hoop and shear stresses, $\sigma_{\theta\theta}(r), \sigma_{r\theta}(r) \approx 0$.

Similarly to Subsec.~{\ref{subsec:II-wrinkling}}, we proceed by considering the displacement and stress fields in the three parts of the sheet: {\emph{(i)}} $R<r<R_{\mathrm{sheet}}$ -- where the sheet is nearly planar, but (unlike Subsec.~{\ref{subsec:II-wrinkling}}) it develops radial wrinkles in $R<r<L_O$ and is axisymmetrically deformed only at $L_O<r<R_{\mathrm{sheet}}$, where both radial and hoop stress components are tensile; {\emph{(ii)}} $L_I<r<R$ -- where the suspended sheet is wrinkled; {\emph{(iii)}} $r<L_I$ -- where the suspended sheet is unwrinkled and both hoop and radial stresses are tensile.
For the last two parts, we notice that the displacement and stress are given by expressions identical to their counterparts in Subsec.~{\ref{subsec:II-wrinkling}}, namely, Eqs.~(\ref{eq:wrinkled-I}) and the BCs~(\ref{eq:BC-sliding-wrinkling-I}) for the nonlinear FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}) in the unwrinkled core, albeit with a different triplet of constants $\Psi(1),\tilde{a}, \tilde{L}_{I}$, that must be determined by matching the radial displacement and stress at the hole's edge
with the wrinkled portion of the sheet at $r>R$. Thus, among the three equations that specify the constants $\Psi(1),\tilde{a}, \tilde{L}_{I}$, the two equations that reflect these continuity conditions
are identical to their counterparts in Eq.~(\ref{eq:twomore}).

In order to find the remaining equation that relates the constants $\Psi(1),\tilde{a}, \tilde{L}_{I}$, we turn to discuss the exterior zone,
$r>R$. Once again, we find a direct mapping
to the Lam\'e problem of an annulus under co-axial, co-planar tensile loads, $\sigma_{rr}(R) = \Psi(1)\cdot\sigma_0$ and $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$. For ${\cal R} \gg 1$ and $\Psi(1) >2$ (for which the Lam\'e solution, Eq.~(\ref{eq:Lame-stress}) is unstable), the TFT solution is given by \cite{Davidovitch11}:
\begin{gather}
{R\!<\!r\!<\!L_O} \left\{
\begin{array}{c}
 \sigma_{rr} (r) \!=\! \sigma_{rr}(R)\frac{R}{r}       \\
    \sigma_{\theta\theta} (r) \!=\! 0
\end{array} \right.
\label{eq:Lame-TFT-2-a}  \\
{L_O\!<\!r\!<\!R_{\mathrm{sheet}}}
\left\{
\begin{array}{c}
 \sigma_{rr} (r) \!=\! \sigma_0 + \left(\sigma_{rr}(L_O) - \sigma_0\right)\frac{L_O^2}{r^2}      \\
  \sigma_{\theta\theta} (r) \!= \!\sigma_0 - \left(\sigma_{rr}(L_O) - \sigma_0\right)\frac{L_O^2}{r^2}
\end{array}
\right.
\label{eq:Lame-TFT-2-b}
\end{gather}
\begin{equation}
{\rm where:}\  L_O = \frac{\sigma_{rr}(R)}{2\sigma_0}R=\frac{\Psi(1)}{2}R  \ ,
\label{eq:Lame-TFT-2-c}
\end{equation}
and
the radial displacement at the wrinkled zone, $R<r<L_O$, is given by:
\begin{equation}
{\rm u}_r(r) = r\frac{\sigma_{rr}(r)}{Y} [- \nu - \log(\frac{L_O}{r})] \ .
\label{eq:dis-stress-Lame-2}
\end{equation}
Comparing Eq.~(\ref{eq:dis-stress-Lame-2}) with its counterpart, Eq.~(\ref{eq:dis-stress-Lame}) in Subsec.~\ref{subsec:II-wrinkling}, reveals a dramatic effect associated with the expansion of wrinkles on the supported part upon increasing indentation depth, $\tilde{\delta}$. While Eq.~(\ref{eq:dis-stress-Lame}) shows that ${\rm u}_{\rm r}(R)$ is proportional to the radial stress at the hole's edge, $\sigma_{rr}(R)$, Eqs.~(\ref{eq:Lame-TFT-2-c},\ref{eq:dis-stress-Lame-2}) show that in the presence of wrinkles
the ratio ${\rm u}_{\rm r}(R)/\sigma_{rr}(R)\propto\log(\sigma_{rr}(R)/\sigma_0)$.
As we will show now, this effect has a strong impact on indentation mechanics, associated with
the continuity equation for radial displacement at the hole's edge:
\begin{equation}
\Psi(1)\cdot \log\frac{\Psi(1)}{2\tilde{L}_{I}} - \frac{1}{2}\tilde{a}^2\cdot(1-\tilde{L}_{I}) = 0 \ ,
\label{eq:thirdmatching-1}
\end{equation}
which follows from Eqs.~(\ref{eq:dis-stress-Lame-2}) and (\ref{eq:ur-core})).

Let us inspect Eq.~(\ref{eq:thirdmatching-1}), contrasting it with its counterpart, Eq.~(\ref{eq:thirdmatching}) in Subsec.~\ref{subsec:II-wrinkling}.
Considering the asymptotic limit $\tilde{\delta} \to \infty$, one may easily notice that a solution of the form $\Psi(1) \sim \tilde{\delta}^2 \ , \ \tilde{a} \sim \tilde{\delta} \ ,  \ \tilde{L}_{I} \sim O(1)$ is consistent with Eq.~(\ref{eq:thirdmatching}), but not with Eq.~(\ref{eq:thirdmatching-1}). Instead, a consistent asymptotic solution of Eq.~(\ref{eq:thirdmatching-1}) has the form: $\Psi(1) \sim \tilde{\delta}^2/\log(\tilde{\delta}) \ , \ \tilde{a} \sim \tilde{\delta} \ , \ \tilde{L}_{I} \sim 1/\log(\tilde{\delta})$. Obtaining the numerical values of the pre-factors in these asymptotic relations requires the use of Eqs.~(\ref{eq:wrinkled-I},\ref{eq:BC-sliding-wrinkling-I},\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}), and the detailed calculation is described in App.~\ref{sub-app:wrinkling}.

The results are shown in the green curves in Figs.~\ref{fig:response-1},\ref{fig:stress-edge} and Figs.~\ref{fig:hoop-stress-profile}-\ref{fig:shape}. We note that the presence of wrinkles on the substrate underlies a sub-cubic asymptotic response, namely $F/\delta^3 \sim 1/\log(\tilde{\delta}) \to 0$ as $\tilde{\delta} \to \infty$, reflecting a logarithmic suppression of the radial stress at the hole's edge with respect to the bare indentation-induced stress: $\sigma_{rr}(R) \!\sim\! \tfrac{1}{\log\tilde{\delta}}\!\cdot\! Y\!\cdot\!(\delta/R)^2$. The invasion of wrinkles into the supported zone of the sheet affects strongly also the displacement field, where the slope at the hole's edge now approaches asymptotically the ``natural'' cone angle: $a \to \tfrac{\delta}{R}\cdot[1- O(1/\log\tilde{\delta})]$, and the size of the unwrinkled core vanishes, $L_I = R\cdot \tilde{L}_{I} \sim R/\log\tilde{\delta}$, as is described in the first row of Table II.

\begin{figure}
\includegraphics[width=0.5\textwidth]{Fig8.pdf}
\caption{A side view of the deformation in the suspended part for $\tilde{\delta} =10$. Radial distances are normalized by the hole's radius $R$, and vertical distances are normalized by the indentation depth $\delta$. The colors represent the same types of BCs as in previous figures.}
\label{fig:shape}
\end{figure}

\subsection{The geometric limit: pseudo-linear response \label{subsec:geo-limit}}}
In the preceding section we saw that if the sheet-substrate attachment is sufficiently weak ($\beta \ll 1$),  radial wrinkles expand in the supported part of the sheet, occupying an annular zone whose external radius, $L_O \sim R \cdot \tilde{\delta}^2 /\log\tilde{\delta}$. If $\tilde{\delta}$ is sufficiently large, wrinkles approach the edge of the sheet, causing yet another dramatic change in the distribution of stress in the sheet and its response to the indentation force. (A similar phenomenon has been found for the indentation of a floating sheet \cite{Vella18b,Ripp20}). For a given value of the parameter ${\cal R}$, our numerical results in the preceding section allow us to estimate the value $\tilde{\delta}^{**}$ at which wrinkles reach the far edge:
\begin{gather}
{\cal R} \approx 0.12 (\tilde{\delta}^{**})^2 \log\tilde{\delta}^{**}  \ \  \Rightarrow \nonumber \\
\tilde{\delta}^{**}({\cal R}) \approx 2.07 \sqrt{{\cal R}\log{\cal R}} \cdot [1 +
O(\log{\cal R})]  \ .
\label{eq:delta-sttd}
\end{gather}
For $\tilde{\delta} > \tilde{\delta}^{**}$, the supported part of the sheet is fully wrinkled, and the stress field for any $R<r<R_{\mathrm{sheet}}$ is described by
Eq.~(\ref{eq:Lame-TFT-2-a}). Together with the BC $\sigma_{rr}(R_{\mathrm{sheet}})=\sigma_0$, we find that:
\begin{equation}
\sigma_{rr}(r) = \sigma_0 \frac{R_{\mathrm{sheet}}}{r} \ \ ; \ \ \sigma_{\theta\theta}(r) = 0 \ ,
\label{eq:TFT-22}
\end{equation}
as we described already in the introductory section~\ref{subsec:Heuristic}. In this regime, the value of the unknown $\Psi(1)$ is directly given by Eq.~(\ref{eq:TFT-22}):
\begin{equation}
\Psi(1) = {\cal R} \ ,
\label{eq:thirdmatching-1a}
\end{equation}
and the deformed state is fully described by solving the FvK equations~(\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}) with the BCs~(\ref{eq:BC-sliding-wrinkling-I}), along with replacing Eq.~(\ref{eq:thirdmatching-1}) by (\ref{eq:thirdmatching-1a}), and the two additional equations in (\ref{eq:twomore}). The solution of these equations is described in App.~\ref{sub-app:wrinkling}. We note that  this solution merely determines the numerical pre-factors in the scaling laws we already found in Subsec.~\ref{subsec:Heuristic}, specifically the pseudo-linear response, $F(\delta) \propto \delta$, Eq.~(\ref{eq:pseudo-linear-1}, with $\gamma_{\rm eff} = \sigma_0$).

The results of this calculation, for a dimensionless indentation depth $\tilde{\delta} = 10$ and ${\cal R} = 6$ (such that $\tilde{\delta} > \tilde{\delta}^{**}({\cal R})$, are shown through the red curves in Figs.~\ref{fig:response-1},\ref{fig:stress-edge} and Figs.~\ref{fig:hoop-stress-profile}-\ref{fig:shape}.
As we noted already in Subsec.~\ref{subsec:Heuristic}, the pseudo-linear response reflects an asymptotically-isometric mechanics, whereby the indentation force
``decouples'' from the stretching modulus of the sheet, transmitting work to the puller at the far edge of the sheet, $r=R_{\mathrm{sheet}}$.
Echoing an observation made already for indenting floating polymer sheets \cite{Vella18b,Ripp20}, our results show that after wrinkles reach the far edge the shrinkage of the tensile core zone with $\tilde{\delta}$ becomes much more pronounced ($\tilde{L}_{I} \sim \tilde{\delta}^{-2}$ {\emph{vs.}}
$\tilde{L}_{I} \sim 1/\log\tilde{\delta}$ for $\tilde{\delta} <\tilde{\delta}^{**}({\cal R})$).
In the asymptotically isometric regime, $\tilde{\delta} \gg \tilde{\delta}^{**}({\cal R})$ the suspended portion
approaches the shape of a perfect cone, with
a slope $\frac{\delta}{R}$, superimposed with radial wrinkles.
\\

\section
The role of sheet-substrate attachment \label{sec:residual_compression}}
In the previous section we avoided the need to address {\emph{explicitly}} the effect of sheet-substrate attachment
by considering the two opposite limits
of strong and weak attachment, namely, the parameter regimes {\emph{(i)}} ($\beta \!\gg\! \tilde{\delta}^2\!\gg\! 1$) and {\emph{(ii)}} ($\beta \!\ll\! 1$), respectively.
Notwithstanding the striking difference between these regimes,
(compare brown {\emph{vs.}} green and red curves in Figs.~\ref{fig:response-1}-\ref{fig:stress-edge} and Figs.~\ref{fig:hoop-stress-profile}-\ref{fig:shape}),
in each of them the mechanical response is not affected by the actual values of the bending modulus $B$ and stiffness $K_{\rm sub}$, but only by the tensile load $\sigma_0$ exerted at the far edge, the stretching modulus $Y$, and the indentation depth $\delta$ (as well as ${\cal R}$).
In contrast, in the intermediate regime {\emph{(iii)}}, $1 \!\ll \!\beta\! \ll \!\tilde{\delta}^2$,
the stress and indentation force depend explicitly on $B$ and $K_{\rm sub}$.

In order to elucidate
this distinction
let us consider a narrow annulus
of radius $r$ as an elastic ring of bending modulus $B$ that is forced to contract due to radial displacement ${\rm u}_{\rm r}(r)<0$ (which is given for each parameter regime by the corresponding expressions in Subsecs.~\ref{subsec:II-wrinkling},\ref{subsec:both-wrinkling},\ref{subsec:geo-limit}).
If the contracted ring is forced to retain a circular shape, it must acquire
a ``bare'' hoop strain, $\varepsilon_{\theta\theta}= \tfrac{{\rm u}_{\rm r}(r)}{r}\!<\!0$, and thereby a compressive stress,
\begin{gather}
\sigma_{\theta\theta}^{(bare)}(r) \approx Y \frac{{\rm u}_{\rm r}(r)}{r} \sim -Y(\delta/R)^2 \ ,
\label{eq:bare-stress}
\end{gather}
and correspondingly an energetic penalty $\sim Y  ({\rm u}_{\rm r}(r)/r)^2$. If out-of-plane deflections are allowed,
the ring may respond as an {\emph{elastica}} -- developing wrinkles of wavelength $\lambda$ and amplitude $A$, such that $(\pi A/\lambda)^2 \approx -\tfrac{{\rm u}_{\rm r}(r)}{r}$.

Such a deformation retains the
arclength nearly intact, suppressing the hoop stress to a residual value \cite{Paulsen16,Davidovitch19}.
\begin{gather}
\sigma_{\theta\theta}^{(res)}(r) \approx - 2B/\lambda^2 \ ,
\label{eq:res-stress}
\end{gather}
whose magnitude will be shown to be much smaller than $\sigma_{\theta\theta}^{(bare)}(r)$.
The
wavelength $\lambda$ and consequently the residual hoop stress, is determined by a ``local $\lambda$ law'' \cite{Cerda03,Paulsen16}:
\begin{gather}
\lambda \approx 2\pi (B/K_{\rm eff})^{1/4}   \ ,
\label{eq:local-lam-law-1}
\end{gather}
where
$K_{\rm eff}$ is an
``effective stiffness'', which may be associated with the resistance of the supporting substrate ($K_{\rm eff} \sim K_{\rm sub}$) or with the presence of radial tension that resists a large wrinkle amplitude ($K_{\rm eff} \sim  \sigma_{rr}(r)/r^2$).
Implementing this rule we find different values of $\lambda$
(and consequently the residual stress and energy) in the suspended and supported parts of the sheet
\begin{eqnarray}
r\!>\!R\! &:&  \ \ \  \lambda \sim \left({B}/{K_{\rm sub}}\right)^{1/4}   \nonumber \\
r<R\! &:&  \ \ \ \lambda \sim (\frac{BR^4}{Y \delta^2})^{1/4}  \sim R\sqrt{t/\delta} \
\label{eq:lambda-1}
\end{eqnarray}
\footnote{The second line of Eq.~(\ref{eq:lambda-1}) is valid for $\tilde{\delta} <\tilde{\delta}^{**}(\beta)$, such that $\sigma_{rr}(R) \sim Y(\delta/R)^2$. For $\tilde{\delta} > \tilde{\delta}^{**}(\beta)$, $\lambda$ at the suspended part is estimated by substituting in Eq.~(\ref{eq:local-lam-law-1}) $K_{\rm eff} \sim \sigma_0 {\cal R}/R^2$.
}.
For the suspended part, $r\!\!<\!\!R$, Eqs.~(\ref{eq:bare-stress},\ref{eq:res-stress},\ref{eq:lambda-1}) show that the residual, wrinkle-induced hoop compression $\sigma_{\theta\theta}^{(res)}                                                                                (r)$
is much smaller than its counterpart $\sigma_{\theta\theta}^{(bare)}(r) \sim Y(\delta/R)^2$, and therefore the formation of wrinkles is energetically favorable in $r<R$, regardless of the value of $\beta$. Turning now to the supported part, and addressing first the parameter regimes {\emph{(i)}} ($\beta \!\gg\! \tilde{\delta}^2 \! \gg\! 1$), and {\emph{(ii)}} ($\tilde{\delta} \!>\! \tilde{\delta}_c  \ \& \ \beta\!\ll\! 1$),
an analogous comparison of $\sigma_{\theta\theta}^{(res)}(r)$ and $\sigma_{\theta\theta}^{(bare)}(r) \sim Y(\delta/R)^2$ yields precisely the same conclusion we reached already in Subsec.~\ref{subsec:threshold-2}, namely, the stress in the supported part is given by the axisymmetric Lam\'e solution in the former regime and by the TFT solution in the latter. However,

when inspecting regime {\emph{(iii)}}, $\tilde{\delta}^2 \!\gg\! \beta \!\gg \!1$, we find that the supported portion, $r\!\!>\!\!R$, consists of a zone close to the hole's edge, where $|\sigma_{\theta\theta}^{(res)}(r)|\!\!\ll\!\! |\sigma_{\theta\theta}^{(bare)}(r)|$, and another zone, away from the hole's edge, where $|\sigma_{\theta\theta}^{(res)}(r)|\!\!\gg\!\! |\sigma_{\theta\theta}^{(bare)}(r)|$. This observation reflects the complexity of the mechanical response in this parameter regime, where
the value of the residual hoop compression, $\sigma_{\theta\theta}^{(res)}(r)$, must be taken {\emph{explicitly}} into account through Eqs.~(\ref{eq:res-stress}, \ref{eq:lambda-1}), in order to reliably evaluate the stress field
and thereby the indentation force.    \\

One may find the stress and indentation force in regime {\emph{(iii)}}
by applying a generalized version of tension field theory \cite{Davidovitch19}. Rather than neglecting the contribution of the residual hoop compression to the radial stress altogether,
Eq.~(\ref{eq:res-stress})
is taken as
a non-homogenous source
in the radial force balance equation~(\ref{eq:FvK-dim-2}), yielding for $r>R$:
\begin{equation}
\psi(r)  = \psi_0 - 2\sqrt{K_{\rm sub} B}{r} \ \ \ \ \Rightarrow  \ \ \ \Psi(\rho) = \Psi_0 - \beta \rho \ ,
\label{eq:Airy-new}
\end{equation}
where $\Psi_0$ is a constant determined through matching conditions with the unwrinkled zones at $r<L_I$ and $r>L_O$, similarly to the analysis in Sec.~\ref{sec:clamping_sliding}. The Airy potential (\ref{eq:Airy-new}) describes a {\emph{bending-induced}} radial tension \cite{hure12,Davidovitch19,Tobasco20},
which can be conveniently expressed as:
\begin{equation}
\sigma_{rr}(r)
=  (\sigma_{rr}(L_O) + 2\sqrt{BK_{\rm sub}})\frac{L_O}{r} - 2\sqrt{K_{\rm sub} B} \ .
\label{eq:srr-new}
\end{equation}

Recalling that we focus here on regime {\emph{(iii)}}, $1\ll \beta \ll \tilde{\delta}^2$, and expecting that $\sigma_{rr}(L_O) \!\sim\! \sigma_0$ (since the region $r\!>\!L_O$ is under nearly isotropic tension),  Eq.~(\ref{eq:srr-new}) can be simplified
in $R<r<L_O$ to $\sigma_{rr}(r) \approx 2\sqrt{BK_{\rm sub}} L_O/r$.
Contrasting this simplified expression with Eq.~(\ref{eq:Lame-TFT-2-a}) or Eq.~(\ref{eq:TFT-22}), we notice that
the stress field and thereby the indentation force in regime {\emph{(iii)}} may be determined in an analogous manner to the analysis of regime {\emph{(ii)}} in Subsecs.~\ref{subsec:both-wrinkling},\ref{subsec:geo-limit}, upon replacing in the definition of the dimensionless variable $\Psi$, Eq.~(\ref{eq:dimensionless-var}):
\begin{equation}
\sigma_0 \to \beta\sigma_0 = 2\sqrt{BK_{\rm sub}} \ .
\label{eq:replace-1}
\end{equation}
Hence, at this level of approximation,
expected to be valid up to corrections of $O(\beta^{-1}) \ll 1$, the mechanics in regime {\emph{(iii)}}  $\tilde{\delta}^2 \!\gg \!\beta\!\gg\! 1$, is described by the mechanics of regime {\emph{(ii)}}  ($\beta \!\ll\! 1 \  \& \ \tilde{\delta} \!\gg\! 1$, with the replacement~(\ref{eq:replace-1}). This observation underlies the last row of Table II.

\section{Hindered sliding: clamping the sheet's edge \label{sec:Expansion}}
In the previous sections, we assumed that the far-edge, $r=R_{\mathrm{sheet}}$, is subjected to a fixed tensile load, $\sigma_0$. Here
we consider another basic boundary condition, which may be of interest to an experimenter, whereby the far edge is clamped. Mathematically, this amounts to replacing the BC $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$, with:
\begin{equation}
{\rm u}_{\rm r}(R_{\mathrm{sheet}}) =(1-\nu)\cdot \frac{\sigma_0}{Y}\cdot R_{\mathrm{sheet}} \ ,
\label{eq:new-BC-Rf}
\end{equation}
where $\sigma_0$ is now understood as an isotropic pre-tension in the sheet {\emph{prior}} to clamping its far edge, $r=R_{\mathrm{sheet}}$ (and prior to indenting its center) \cite{Vella18}. Clamping the sheet at its far edge hinders its sliding inwards,
which is necessary to release the radial strain induced by indentation. Thus, for a given $\tilde{\delta} \gg 1$, the in-plane stress in this version of the problem
is larger in comparison to a sheet under fixed tensile load, $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$, and so is the indentation force.
This effect is elucidated by contrasting the corresponding versions of the
Lam\'e problem.
In the first version, which was the basis for our analysis in the preceding sections, the far edge is under a given radial tension, $\sigma_{rr}(R_{\mathrm{sheet}}) = \sigma_0$, but otherwise is free to slide on the substrate (${\rm u}_{\rm r}(R_{\mathrm{sheet}})< 0$),
the stress field of the planar, unwrinkled state, is given by Eq.~(\ref{eq:Lame-stress}), and the tension field solution of the wrinkled state is given by Eqs.~(\ref{eq:Lame-TFT-2-a}-\ref{eq:Lame-TFT-2-c}). In the second version of the Lam\'e problem, the BC at the far edge is given by Eq.~(\ref{eq:new-BC-Rf}), the stress field of the planar (unwrinkled) state is:
\begin{widetext}
\begin{eqnarray}
\sigma_{rr}(r) &=& \frac{{\cal R}^2}{{\cal R}^2 (1-\nu)
 + (1+\nu)} \cdot \left\{
[\sigma_0 (1-\nu) + \sigma_{rr}(R)(1+\nu)]
+ [(1-\nu)[\sigma_{rr}(R) - \sigma_0](\frac{R}{r})^2 \right\}
\label{eq:Lame-srr-clamped} \\
&\approx& \frac{1}{1-\nu}
 \cdot \left\{
[\sigma_0 (1-\nu) + \sigma_{rr}(R)(1+\nu)]
+ [(1-\nu)[\sigma_{rr}(R) - \sigma_0](\frac{R}{r})^2 \right\}
\nonumber \\
\sigma_{\theta\theta} (r) &=& \frac{{\cal R}^2}{{\cal R}^2 (1-\nu) \!+\! (1+\nu)} \cdot \left\{
[\sigma_0 (1-\nu) \!+\! \sigma_{rr}(R)(1+\nu)]
\!-\! [(1-\nu)[\sigma_{rr}(R) \!-\! \sigma_0](\frac{R}{r})^2 \right\}  \ .
\label{eq:Lame-sqq-clamped} \\
&\approx& \frac{1}{1-\nu}
\cdot \left\{
[\sigma_0 (1-\nu) \!+\! \sigma_{rr}(R)(1+\nu)]
\!-\! [(1-\nu)[\sigma_{rr}(R) \!-\! \sigma_0](\frac{R}{r})^2 \right\}  \ ,  \nonumber
\end{eqnarray}
\end{widetext}
where the second lines in the above equations are valid for ${\cal R} \gg 1$.
As a result, the TFT solution is characterized by a compression-free stress in the wrinkled zone:
\begin{equation}
R<r<L_O : \left\{
\begin{array}{c}
   \sigma_{rr} (r) = \sigma_{rr}(R)\frac{R}{r}    \\
    \sigma_{\theta\theta} (r) = 0
\end{array}
\right.
\label{eq:Lame-TFT-2-clamped}
\end{equation}
with
\begin{equation}
\resizebox{0.4\textwidth}{!}{$
L_O \!=\!R \left(
\sqrt{\frac{1-\nu}{1+\nu}} \!\cdot\! \sqrt{1\!+\! (\frac{{\cal R}}{\Psi(1)})^2\frac{1-\nu}{1+\nu}}\!\cdot\! {\cal R}  -
\frac{{\cal R}^2}{\Psi(1)}\frac{1-\nu}{1+\nu} \right)  $}
\label{eq:Lame-TFT-2-clamped-Lo}
\end{equation}
\vspace{-1.0cm}
$$
\approx \frac{1}{2}R \Psi(1)    \ \ \ ({\rm for} \ \ {\cal R} \gg 1)  \ , $$
where $\Psi(1) = \sigma_{rr}(R)/\sigma_0$, and the stress components in the unwrinkled zone, $L_O<r<R_{\mathrm{sheet}}$, are given by
Eqs.~(\ref{eq:Lame-srr-clamped},\ref{eq:Lame-sqq-clamped}) upon replacing:
$R \!\to\! L_O, \sigma_{rr}(R) \!\to\! \sigma_{rr}(L_O) \!= \!\sigma_{rr}(R)\!\cdot\! R/L_O$, and ${\cal R} \!\to\! R_{\mathrm{sheet}}/L_O \!=\! {\cal R}\!\cdot\! R/L_O$.

The primary effect of the BC~(\ref{eq:new-BC-Rf}) is elucidated by considering a fixed ${\cal R}\!\gg\! 1$, and using the above expressions to evaluate $\sigma_{rr}(R_{\mathrm{sheet}})$ for $\Psi(1) \!= \!\sigma_{rr}(R)/\sigma_0 \!\to\! \infty$.
For both planar state and wrinkled state, we find that the far-edge stress $\sigma_{rr}(R_{\mathrm{sheet}})$ it proportional to the stress at the hole's edge $\sigma_{rr}(R)$. More specifically, we find that for the planar state $\sigma_{rr}(R_{\mathrm{sheet}})/\sigma_{rr}(R) \propto {\cal R}^{-2}$,
whereas for the wrinkled state $\sigma_{rr}(R_{\mathrm{sheet}})/\sigma_{rr}(R) \propto {\cal R}^{-1} $.
This means that in order to keep the far edge from sliding inwards under the influence of the large radial stress $\sigma_{rr}(R)$ that pulls at the inner edge, the clamp must exert a comparable radial load on the far edge, hence $\sigma_{rr}(R_{\mathrm{sheet}}) \sim \sigma_{rr}(R) \gg \sigma_0$.
This observation is rather intuitive, indicating
that the elastic energy needed to deform a sheet clamped at its far edge
is much larger than the energy required to deform a sheet whose far edge is free to slide. As a consequence, the indentation force $F(\delta)$
is larger in comparison to the response we found in the preceding sections for a sheet subjected to a fixed boundary load.

We find the indentation force $F(\delta)$ by following the tracks of our analysis in Subsecs.~\ref{subsec:II-wrinkling},\ref{subsec:both-wrinkling},
assuming the sheet is wrinkled in an azimuthally-confined zone, $L_I<r<L_O$ (with $L_I <R$ and $L_O >R$) and unwrinkled in $r<L_I$ and $L_O<r<R_{\mathrm{sheet}}$.
Considering the first two zones, we notice that the displacement and stress are given by expressions identical to Eqs.~(\ref{eq:wrinkled-I}) and the BCs~(\ref{eq:BC-sliding-wrinkling-I}) for the nonlinear FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}) in the unwrinkled core, albeit with a different triplet of constants $\Psi(1),\tilde{a}, \tilde{L}_{I}$, that must be determined by matching the radial displacement at the hole's edge ${\rm u}_{\rm r}(R)$ with the wrinkled state at the exterior of the hole. Thus, exactly as we found
in Subsec.~\ref{subsec:both-wrinkling}, two equations among the three
that specify the constants $\Psi(1),\tilde{a}, \tilde{L}_{I}$,
are identical to their counterparts in Eq.~(\ref{eq:twomore}), and the third equation reflects a continuity of the radial displacement at $r=R$. Employing Eqs.~(\ref{eq:Lame-srr-clamped}-\ref{eq:Lame-TFT-2-clamped-Lo}) and Eq.~(\ref{eq:ur-core}) we obtain an equation that replaces Eq.~(\ref{eq:thirdmatching-1}):
\begin{widetext}
\begin{equation}
\Psi(1)\cdot \log\left(-\frac{{\cal R}^2}{\Psi(1)\tilde{L}_{I}}\frac{1-\nu}{1+\nu} \ + \frac{{\cal R}}{\tilde{L}_{I}}\sqrt{\frac{1-\nu}{1+\nu}}\sqrt{1+\frac{1-\nu}{1+\nu}(\frac{{\cal R}}{\Psi(1)})^2} \right)
- \frac{1}{2}\tilde{a}^2\cdot(1-\tilde{L}_{I}) = 0 \ ,
\label{eq:thirdmatching-2}
\end{equation}
\vspace{-0.5 cm}
$$ \Longrightarrow \Psi(1)\cdot \log\left( \frac{{\cal R} \Psi(1)}{2 \tilde{L}_{I}} \right) - \frac{1}{2}\tilde{a}^2\cdot(1-\tilde{L}_{I}) \approx 0  \ \ \
({\rm for} \ {\cal R} \gg 1)
\ .  $$
\end{widetext}

The values of the unknowns $\Psi(1), \tilde{L}_{I}, \tilde{a}$, for any given $\tilde{\delta} > 3.3$ are obtained from the numerical solution of the three algebraic equations that are derived from Eqs.~(\ref{eq:twomore},\ref{eq:thirdmatching-2}) and the exact solutions of the FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}), under the BCs~(\ref{eq:BC-sliding-wrinkling-I}); see details in App.~\ref{sub-app:wrinkling}.
As anticipated by the above discussion, we notice that if the far-edge is clamped, the asymptotic response at large indentation depth, $\tilde{\delta} \to \infty$, is $F \sim (Y/R^2)\cdot \delta^3$, hence the system does not reach the extreme wrinkle-assisted softening obtained upon exerting a fixed load at the far edge. Nevertheless, the asymptotic value of the constant, $F/\delta^3$, scales as $1/\log({\cal R})$ (second row of Table II), so that as ${\cal R}$ is increased, the wrinkle-induced suppression of the indentation force becomes more and more effective.

\section
Wrinkling {\emph{vs.}} delamination} \label{sec:wrink-delam}
In our model we assume that relieving compression in the supported part of the sheet does not require the formation of delaminated zones, in which the sheet-substrate distance $d$ exceeds the width of the VdW potential well (Fig.~\ref{fig:schem-slide}d), but merely tiny deviations of $d$ from the thermodynamic equilibrium value $d_{\rm min}$. The crucial distinction between these deformation types stems from the
respective
energetic costs (per area) of sheet-substrate attachment:
\begin{gather}
{\rm delamination}:  u_{\rm sub} \!\approx\! V({d_{\rm min}})
\label{eq:VdWmodel-1} \\
{\rm Zhang\!-\!Witten}:
u_{\rm sub}(d) \!\approx\! \tfrac{1}{2}V''({d_{\rm min}})\!\cdot\! (d\!-\!d_{\rm min})^2  \label{eq:VdWmodel-2}
\end{gather}
With the Zhang-Witten stiffness, $K_{\rm sub} = V''({d_{\rm min}})$,
a rigid substrate that supports a thin sheet is merely an example of a ``Winkler foundation'' \cite{TimoshenkoBook}, hence the response of the sheet to compression is analogous to other examples of this basic model, such as a sheet floating on a liquid bath (where $K_{\rm sub} = \rho_{liq} g$ with $\rho_{liq}$ being the liquid's mass density). For Winkler-like problems, planar deformations are unstable to wrinkling -- periodic undulations characterized by a {\emph{single wavelength}} $\lambda \sim (B/K_{\rm sub})^{1/4}$ (see Eq.~\ref{eq:local-lam-law-1}) -- which emerges through a supercritical (second order) instability of the planar state, not involving any energy barrier.

In contrast, the finite, $d$-independent energy $V(d_{\rm min})$ associated with delamination, Eq.~(\ref{eq:VdWmodel-1}), which one may view as a surface energy penalty, entails a strictly different instability of the planar state. This instability is sub-critical (first order), and therefore requires the crossing of an energy barrier, which in turn gives rise hysteresis loops. Furthermore, the basic deformation mode \cite{Wagner12,Davidovitch20} is a {\emph{single delaminated zone}}, which may accommodate any excess length by increasing the sheet-substrate distance $d$ without further energy cost, as is indicated by Eq.~(\ref{eq:VdWmodel-1}), rather than by forming multiple delamination zones. Even though periodic delamination patterns have been observed under certain circumstances (such as uniaxial compression of a sheet attached to compliant substrate \cite{Vella09}), those patterns are characterized by two length scales, whereby the width of each delaminated zone is much smaller than the distance between them (where the sheet remains fully laminated). Hence, even if the  indentation-induced hoop compression leads to delamination instability,
the number of blisters at a given distance $r$ should be $\ll 2\pi r/\lambda$, where $\lambda$ is the average width of an individual blister. This suggests that a recent attempt to describe such a delamination pattern by a wrinkling-like sinusoidal profile, characterized by single wavelength $\lambda$ \cite{Dai2020}, is nonphysical.

In order to determine which of the two deformation types, described by Eqs.~(\ref{eq:VdWmodel-1}) and (\ref{eq:VdWmodel-2}), is likely to relieve hoop compression in a given indentation experiment,
we note two necessary conditions for a wrinkle pattern to be physically realizable.

{\emph{(i)}} The wrinkle wavelength $\lambda$ (Eq.~\ref{eq:local-lam-law-1}, with $K_{\rm eff}=K_{\rm sub}$, Eq.~\ref{eq:KsubWitten}) must exceed the length $\ell_{bend} = \sqrt{B/Y}$, otherwise the bending energy would be too large, rendering wrinkles energetically unfavorable. In terms of the parameters of our model this condition reads:
\begin{equation}
\ell_{bend} \ll \ell_{VdW}  \ ,
\label{eq:cond-wrink-1}
\end{equation}
where we defined the length scale:
\begin{equation}
\ell_{VdW} \equiv \sqrt{Y/V''(d_{min})} \ ,
    \label{eq:ell-vdw}
\end{equation}

{\emph{(ii)}} The wrinkle amplitude $d$ must not exceed a length $d_{max}$ above which the sheet ``escapes'' from the attractive zone of the VdW potential (see schematic Fig.~\ref{fig:schem-slide}d), and the energetic cost transitions from Eq.~(\ref{eq:VdWmodel-2}) to Eq.~(\ref{eq:VdWmodel-1}). Noting that the ratio between the wrinkle amplitude and wavelength is ``slaved'' to the excess hoop length,
$(|d-d_{min}|/\lambda)^2 \sim -{\rm u}_{\rm r}/r$ (such that the wrinkly undulations ``waste'' just the right arclength necessary to suppress hoop compression \cite{Davidovitch11}), and using the estimate ${\rm u}_{\rm r} \sim - \delta^2/R$ (Eq.~\ref{eq:radial-scale}), we obtain the second condition:
\begin{equation}
\frac{\delta}{R} \ll \frac{|d_{max}-d_{min}|}{\sqrt{\ell_{bend} \cdot \ell_{VdW}}} \ .
    \label{eq:cond-wrink-2}
\end{equation}
The two conditions~(\ref{eq:cond-wrink-1},\ref{eq:cond-wrink-2}) define a parameter regime in which we expect the wrinkle patterns assumed in our model to be a feasible, energetically-favorable mechanism for relaxing the hoop compression induced by indentation and sliding. If condition~(\ref{eq:cond-wrink-1}) is violated, an axisymmetric (unwrinkled) deformation in the supported portion of the sheet (SubSec.~\ref{subsec:II-sliding}) is stable against wrinkling, and delamination may occur, through a sub-critical instability, at some large  indentation depth directly from the planar state. If condition~(\ref{eq:cond-wrink-1}) is satisfied, the supported portion of the sheet becomes unstable to wrinkling at $\tilde{\delta}^{*}(\beta)$ (Eq.~\ref{eq:threshold-stst}), and delamination is expected to occur when the indentation depth $\delta$ reaches ${ R (d_{max}-d_{min})}/{\sqrt{\ell_{bend} \cdot \ell_{VdW}}}$.

A crude estimate of the various lengths in the conditions~(\ref{eq:cond-wrink-1},\ref{eq:cond-wrink-2}) may be obtained by assuming
$V''(d_{min}) \sim V(d_{min})/d_{min}^2$, and $ 0.1 {\rm nm} <d_{min} < d_{max} < 1 {\rm nm}$. For Graphene (on SI or BN), we use the values $Y \sim 300 N/m \ , \ \ell_{bend} \sim 0.1 {\rm nm}$, and $V(d_{min}) \sim 0.1 N/m$. With these values,
we find that both conditions~(\ref{eq:cond-wrink-1}) and (\ref{eq:cond-wrink-2}) are satisfied for $\delta/R < 0.1$,
suggesting the relevance of a wrinkle-assisted compression-relieving mechanism for experiments, at least at indentation depths $\delta \lesssim 100$ nm.

\section{Discussion \label{sec:discussion}}

\subsection{The non-perturbative macroscale effect of bending rigidity}
Employing standard TFT (Secs.~\ref{sec:clamping_sliding},\ref{sec:Expansion}) or its recently generalized version (Sec.~\ref{sec:residual_compression}) we showed that, as long as there is
compressive stress somewhere within the indented sheet, the ability to relax it by energetically-inexpensive wrinkles acts to suppress considerably the elastic energy. Our results, summarized in Tables I and II, show that the wrinkle-assisted reduction of elastic energy and the consequent suppression of the indentation force $F(\delta)$ is a {\emph{non-perturbative}} phenomenon, which is not sensitive to the specific value of the bending modulus, but rather stems from its mere smallness ({\emph{i.e. $\epsilon \ll 1$}}). That is, for specific BCs ({\emph{e.g.}} sliding at $r=R$ and a constant tensile load at $r=R_{\mathrm{sheet}}$), we find that the error incurred by ignoring the effect of wrinkles on the indentation force is $O(Y/R^2)\delta^3$, as one can see by comparing the second row of Table I (which ignores the effects of wrinkles, describing a mechanically-unstable state
for $\tilde{\delta} \!>\!\tilde{\delta}_c$) with the third row of Table I or the first two rows of Table II.

While we focused our study on the pointwise indentation problem, the above lesson is general and applies to any situation in which a confining geometry or loading conditions give rise to compressive stress within a thin, highly bendable sheet. One example, which has attracted some interest lately, is the strain induced in a 2D solid sheet, supported on a smooth substrate, by
high-pressure ``bubbles'' confined between the sheet and the substrate \cite{Khestanova16}. Such bubbles cause radial stretching of the sheet around the bubble axis, and -- similarly to the indentation problem (with sliding BCs) -- a hoop compression emerges in the sheet at the vicinity of the bubble's edge. While a wrinkle-assisted suppression of hoop compression may not have a pronounced effect on the bubble's shape or the pressure within it \cite{Khestanova16}, the strain components in the sheet are strongly affected by the presence of wrinkles. This effect, however, has been overlooked in a recent paper \cite{Dai18}, where the authors computed the strain tensor by assuming a mechanically-unstable (unwrinkled, axisymmetric) deformation of the sheet.

\subsection{Beyond ideal mechanics -- substrate roughness and thermal fluctuations}
Our model assumes a smooth, homogeneous substrate, such that the only energetic cost of sliding stems from the consequent hoop compression. From a pure mechanical perspective, a roughness of the substrate may give rise to localized or extended zones in which the sheet is pinned to the substrate, hindering its sliding inwards. A simple, effective-medium-theory approach to incorporate surface roughness into our model may be to replace the {\emph{control parameters}} $\gamma_{\rm eff}$ and ${\cal R}$ in the last two rows of Table II with {\emph{effective parameters}} that account for the excess radial tensile and clamping (away from the hole's edge), associated with the hindrance of sliding. A more thorough study of the effects of surface roughness, as well as thermal fluctuations, on the indentation force, should account for the anomalous elasticity\cite{NP87,AL88,LR92,Radzihovsky18} that has been predicted for 2D solid membranes such as Graphene at room temperature \cite{Blees14,Kosmrlj17,KatsnelsonBook}.

\subsection{Summary}

The main purpose of the
ideal model we introduced in this paper is to elucidate the crucial assumptions one has to make in order to extract the stretching modulus of a suspended sheet from indentation experiments. In this context, the central outcome of our analysis is that
sliding and wrinkling of the sheet affect significantly the commonly-assumed cubic dependence of the indentation force, $F/\delta^3 \propto (Y/R^2)$; the assumption of clamping at the edges of the suspended sheet gives a lower bound to the value of the Young modulus. If the membrane can slide over the non suspended zone, the force required to achieve a given deformation can be significantly lower than in the case of clamping.
This message is illustrated most conspicuously
in the geometry-dominated nature of the pseudo-linear response, Eq.~(\ref{eq:pseudo-linear-1}), where $F/ \delta$ may depend on a pre-tension $\sigma_0$ or a bending-induced tension $\gamma_{\rm eff} = 2\sqrt{BK_{\rm sub}}$, as well as on the radii $R$ (of the hole) and $R_{\mathrm{sheet}}$ (of the whole sheet), but not on the stretching modulus $Y$! Such a stretching-independent
response may be avoided if the attachment to the substrate is sufficiently strong, or if the sheet is clamped at the far edge ($r = R_{\mathrm{sheet}} \gg R$). But also in such cases
sliding and wrinkling have a significant effect on the indentation force, which must be considered
in order to properly extract the stretching modulus $Y$ from the measured response.

Our theoretical model is quite elementary and does not include effects which may be important for experimental set-ups of 2D membranes, such as pinning, spatial disorder, and thermal fluctuations. We suspect that further theoretical progress is required, possibly along the directions outlined above, in order to render our model applicable for a quantitative description of actual experiments. Nevertheless, some basic predictions may be sufficient to test the relevance (or lack thereof) of sliding and wrinkling. Specifically, measuring the slope ($\approx \theta$) of the suspended sheet in the vicinity of the hole's edge may provide a robust, indirect probe for this purpose. A slope that is close to $63\%$ should indicate that the sheet is practically clamped at the hole's edge. A larger slope should indicate a substantial sliding and wrinkling of the sheet in the suspended part and possibly also on the substrate.

Beyond its relevance to metrology
and to studying sliding and wrinkling phenomena, our model highlights the complexity that is often ignored by one's perception of 2D solid membranes as being ``nearly inextensible, highly bendable'' objects, whose resistance to bending can be ignored in analyzing macro-scale, tension-dominated deformations. Instead, our study
illuminates the subtle role played by both stretching and bending rigidity in the response to such external stimuli.

\begin{acknowledgements}
We thank A. Geim, M. Katsnelson, K. Novoselov, D. Vella, and participants of the program ``Geometry and Elasticity of 2D Soft Matter'' at the Kavli Institute for Theoretical Physics Santa Barbara 2016, where we began working on this manuscript,
for many useful discussions. We thank D. Vella for a thorough, critical review of the manuscript.
We acknowledge support by the National Science Foundation under grants NSF-DMR-CAREER-1151780 and NSF-DMR-1822439 (BD), and by the European Commission, under the Graphene Flagship, Core 3, grant no. 881603, and by the grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat and SEV-2016-0686, (Ministerio de Ciencia e Innovacion, Spain) (FG).
\end{acknowledgements}

\begin{appendix}
\section{Boundary conditions at the hole's edge and the negligibly of radial curvature} \label{app:krr-BC}
In order to elucidate the neglect of the radial bending force $B\partial^4z/\partial r^4$ in Eq.~(\ref{eq:FvK-dim-1}), often referred to as a ``membrane approximation'', as well as the BCs at the hole's edge, let us recall that for a sheet with finite (albeit small) bending modulus $B$ the tangent $\hat{t}$ to the sheet's plane must be a continuous function of the radial distance $r$. A discontinuity of $\hat{t}$ implies a divergence of the radial curvature, $\kappa_{rr} \approx |\partial\hat{t}/\partial r|$, and hence an infinite bending energy, regardless of how small $B$ is. In fact, the vicinity of the hole's edge, where the tangent $\hat{t}$ varies sharply, is the only zone where the radial curvature has to be considered, since it is required to regularize this divergence.
Specifically, the characteristic length over which occurs the necessary change from $\hat{t}(r\to R^-) \!=\! \cos\theta \hat{r} \!-\! \sin\theta \hat{z}$ to $\hat{t}(r\to R^+) \!=\! \hat{r}$, is the ``local bendo-capillary'' length \cite{Vella18}:
\begin{equation} \ell_{bc}^* \approx \sqrt{B/\sigma_{rr}(R)}  \ . \label{eq:bendo-cap}
\end{equation}
On one hand, we have that $\ell_{bc}^*\ll R$, since the sheet is highly bendable ({\emph{i.e.}} $\epsilon\ll 1$, see Eq.~(\ref{eq:FvK-nondim-2}) and the following paragraphs); on the other hand we assume $\ell_{bc}^*$ is much larger than the atomic scale (over which the corner in the substrate is ``smoothed out''), see schematic Fig.~\ref{fig:schem-slide}b-c.

In our analysis of the FvK equations, either of the unwrinkled state in Subsec.~\ref{subsec:II-sliding} or the wrinkled state in Subsec.~\ref{subsec:II-wrinkling} and the rest of the paper, we exploited the fact that $\ell_{bc}^* \ll R$, and considered the narrow annulus,
$R-\ell_{bc}^* < r< R$, as a ``boundary layer'', whose energetic cost may be ignored. More precisely, this excess energy can be estimated as
$\sim B(\theta/\ell_{bc}^*)^2 R \sim \sqrt{B Y} \delta^3/R^2$, and an inspection of Tables I-II reveals that it is smaller by a factor $\sqrt{\epsilon}$, Eq.(~\ref{eq:DG-0}), than
the elastic energy evaluated in Secs.~\ref{sec:clamping_sliding}-\ref{sec:Expansion}.
Hence, neglecting the explicit energetic cost of that boundary
layer amounts to evaluating the leading order of the elastic energy (and the indentation force derived from it) in an expansion whose small parameter is $\sqrt{\epsilon}$.
Mathematically, since the radial bending force, $B\partial^4z/\partial^4r$, is significant only in this narrow zone, our analysis has been greatly simplified by omitting this term from the $1^{st}$ FvK Eq.~(\ref{eq:FvK-dim-1}), rendering it -- along with Eq.~(\ref{eq:FvK-dim-2}) -- a coupled set of $2^{nd}$ order ODEs for $\psi(r)$ and $z(r)$, and allowing for a discontinuity of $z'(r)$ at $r=R$.

The boundary layer approach
implies that the radial and vertical components of the displacement may be considered continuous at $r=R$ yielding the BCs~(\ref{eq:matchingur})  and~(\ref{eq:BC-nondim-sliding-axi}iii), while the derivative of the latter is allowed to be discontinuous ($[z'(r)]_{R^-}^{R+} \approx \theta$). At the same time, the mere existence of the boundary layer underlies the continuity of the radial stress component (even though one may naively view it as a violating a force balance in the horizontal direction at $r=R$), as is illustrated in the schematic Fig.~\ref{fig:schem-slide}. We note that these continuity BCs remain valid even if a small portion of the sheet slides vertically in order to gain some surface energy by contacting the hole's walls (contrast panels b and c in Fig.~\ref{fig:schem-slide}), as long the sheet does not get pinned to the substrate. A detailed discussion of this effect will be discussed elsewhere.

\section{General analysis of the unwrinkled core \label{app:axisymmetric}}
Here we describe the steps underlying an analytic solution for an axisymmetric (unwrinkled) solution the nonlinear FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}). This solution, with distinct types of BCs, is used to characterize a purely tensile ``core'' around the indenter, which exists under all various conditions (clamping/sliding at the hole's edge, and various parameter regimes, Eqs.~(\ref{eq:regime-i},\ref{eq:regime-ii})). Our exposition follows closely Ref. \cite{Vella17} and the Supplementary Information of Ref. \cite{Vella15}
\\

We start by integrating the $1^{st}$ FvK equation (\ref{eq:FvK-nondim-1}), and obtain:
\begin{equation}
\Psi \frac{d \zeta}{d \rho} = {\cal F}
\label{eq:integrate-1}
\end{equation}
Next, we introduce the variable transformation \cite{Bhatia68}:
\begin{equation}
\Phi = \rho \Psi \ \ , \ \ \eta = \rho^2  \ ,
\label{eq:transform-1}
\end{equation}
such that: $\Psi = \tfrac{\Phi}{\sqrt{\eta}}$, and $\tfrac{d \Psi}{d \rho} = \left( 2 \tfrac{d \Phi}{d \eta} - \frac{\Phi}{\eta} \right)$.
With this transformation, the $2^{nd}$ FvK equation (\ref{eq:FvK-nondim-2}) becomes:
\begin{equation}
\Phi'' = - \frac{{\cal F}^2}{8 \Phi^2} \ ,
\label{eq:FvK-nondim-tran-2}
\end{equation}
which can be integrated once, obtaining:
\begin{align}
\Phi' &= \frac{\cal F}{2} \frac{\sqrt{1 + A \Phi}}{\sqrt{\Phi}}
\label{eq:FvK-nondim-tran-integrate-2}
\end{align}
where $A$ is a constant of integration. Evaluating Eq.~(\ref{eq:FvK-nondim-tran-integrate-2}) at $\eta=1$, we obtain a first equation that involves the unknowns $\Phi'(1), \Phi(1), A$, and ${\cal F}$:
\begin{align}
\Phi'(1) &= \frac{\cal F}{2} \frac{\sqrt{1 + A \Phi(1)}}{\sqrt{\Phi(1)}}
\label{eq:FvK-nondim-tran-integrate-2-edge}
\end{align}

Integrating now Eq.~(\ref{eq:FvK-nondim-tran-integrate-2}), we obtain an explicit expression between the variable $\eta$ and the function $\Phi(\eta)$:
\begin{align}
\frac{\sqrt{\Phi(1 + A \Phi)}}{A}   - \sqrt{\frac{1+A\Phi}{A^3}} \sinh^{-1}[\sqrt{A\Phi}] &=  \frac{\cal F}{2} \eta  \ ,
\label{eq:FvK-nondim-tran-integrate-22}
\end{align}
(where we used the BC $\Phi(0) = 0$, which is valid for all cases addressed here). Evaluating the above equation at the hole's edge ($\eta = 1$), we obtain a second equation that involves the unknowns $\Phi(1),A$, and ${\cal F}$:
 \begin{align}
\frac{\sqrt{\Phi(1)(1 + A \Phi(1))}}{A}   - \sqrt{\frac{1+A\Phi(1)}{A^3}} \sinh^{-1}[\sqrt{A\Phi(1)}] &=  \frac{\cal F}{2}  \  .
\label{eq:FvK-nondim-tran-integrate-22-edge}
\end{align}

Turning to the integrated form of the $1^{st}$ FvK Eq.~(\ref{eq:integrate-1}), we re-parametrize the function $\zeta(\rho) \to \zeta[\Phi(\eta)]$. With the aid of Eq.~(\ref{eq:FvK-nondim-tran-integrate-2}), and integration (over $\Phi$), we obtain an explicit form for the shape:
\begin{equation}
\zeta(\Phi) - \zeta(0) = \frac{2}{\sqrt{A}} \sinh^{-1}[\sqrt{A\Phi}]  \ .
\label{eq:FvK-nondim-tran-integrate-1}
\end{equation}
\\

Equations~(\ref{eq:FvK-nondim-tran-integrate-2-edge},\ref{eq:FvK-nondim-tran-integrate-22-edge}) constitute  two equations for the four unknowns: $\Phi'(1), \Phi(1), A$, and ${\cal F}$. These two equations are common to all cases we study in this paper. The other two equations must come from the BCs that reflect the various physical conditions discussed in our paper (clamping/sliding at the hole's edge, absence/presence of wrinkles).

Once the four constants ($\Phi'(1), \Phi(1), A ,{\cal F}$) are determined, Eqs.~(\ref{eq:FvK-nondim-tran-integrate-22},\ref{eq:FvK-nondim-tran-integrate-1}) provide explicit expressions for the functions $\Phi(\eta), \zeta(\Phi)$, which can be directly transformed (through Eqs.\ref{eq:transform-1},\ref{eq:define-Airy}) to the shape, $\zeta(\rho)$, and the stress components: $\sigma_{rr}(r),\sigma_{\theta\theta}(\rho)$.

\section{Clamping at the hole's edge}
\label{sub-app:clamping}
For the clamped case, Subsec.~\ref{subsec:II-clamping}, the BCs (\ref{eq:BC-nondim-clamping}) become:
\begin{align}
\eta =0 \ &:  \, \, \, \,  (i ) \, \zeta = -\tilde{\delta} \, \, \, \, \, \, \, \, \,(ii ) \Phi = 0 \nonumber \\
\eta =1 \ &:  \, \, \, \,  (iii) \, \zeta = 0 \, \, \, \, \, \, \, \, \, (iv) \, 2\Phi'  = (1-\nu) + (1+\nu) \Phi  \ .
\label{eq:BC-nondim-trans-clamping}
\end{align}
Among these BCs, {\emph{(ii)}} was used already to obtain Eq.~(\ref{eq:FvK-nondim-tran-integrate-22}). Since the FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}) are invariant under: $\zeta \to \zeta + c$, only the difference $\zeta(1) - \zeta(0)$ can affect the physics, and hence the three remaining BCs in (\ref{eq:BC-nondim-trans-clamping}) give rise to two equations that involve the unknowns ($\Phi'(1), \Phi(1), A,{\cal F}$). The first equation is simply BC {\emph{(iv)}}:
\begin{equation}
2\Phi'(1)  = (1-\nu) + (1+\nu) \Phi(1)  \ ,
\label{eq:BC-iv-clamping}
\end{equation}
and the second equation is obtained by evaluating Eq.~(\ref{eq:FvK-nondim-tran-integrate-1}) at $\Phi(1)$, and
substituting for the difference: $\zeta(\Phi(1))-\zeta(\Phi(0)) = \tilde{\delta}$:
\begin{equation}
\tilde{\delta} = \frac{2}{\sqrt{A}} \sinh^{-1}[\sqrt{A\Phi(1)}]  \ .
\label{eq:FvK-nondim-tran-integrate-1-edge}
\end{equation}

Solving the four algebraic equations (\ref{eq:FvK-nondim-tran-integrate-2-edge},\ref{eq:FvK-nondim-tran-integrate-22-edge},\ref{eq:BC-iv-clamping},\ref{eq:FvK-nondim-tran-integrate-1-edge}), is straightforward ({\emph{e.g.}} using Mathematica's ``FindRoot"), and allows us to obtain the constants
$\Phi'(1), \Phi(1), A,{\cal F}$, as a function of the single dimensionless parameter $\tilde{\delta}$. The response function ${\cal F} (\tilde{\delta})$, the deformed shape, and the stress profile (which are evaluated with the aid of
Eqs.~(\ref{eq:FvK-nondim-tran-integrate-22},\ref{eq:FvK-nondim-tran-integrate-1},\ref{eq:transform-1},\ref{eq:define-Airy}), are shown in the gray curves in Figs.~\ref{fig:response-1}-\ref{fig:shape}.

\section{Sliding (no wrinkling)}
\label{sub-app:sliding}

The BCs that corresponds to an axisymmetric (unwrinkled) state, for which the sheet can slide on the substrate were derived in Subsec.~\ref{subsec:II-sliding}. The difference between clamped-edge and sliding boils down to replacing the BC {\emph{(iv)}} in  Eq.~(\ref{eq:BC-nondim-clamping}) with the corresponding BC in Eq.~(\ref{eq:BC-nondim-sliding-axi}). Hence, the algebraic equations for the four unknowns $\Phi'(1), \Phi(1), A,{\cal F}$ are Eqs.~(\ref{eq:FvK-nondim-tran-integrate-2-edge},\ref{eq:FvK-nondim-tran-integrate-22-edge},\ref{eq:FvK-nondim-tran-integrate-1-edge}), and:
 \begin{equation}
\Phi'(1)  = 1  \ .
\label{eq:BC-iv-sliding-axi}
\end{equation}
The response function ${\cal F} (\tilde{\delta})$, the deformed shape, and the stress profile
that correspond to this solution,
are shown in the blue curves in Figs.~\ref{fig:response-1}-\ref{fig:shape}.

\section{Sliding and Wrinkling
\label{sub-app:wrinkling}}
If the sheet can slide at the hole's edge, hoop compression evolves around the hole's edge and the compressed zone expands upon increasing indentation depth,  in a manner that depends on the sheet-substrate attachment (through the parameter $\beta$, Subsecs.~\ref{subsec:II-wrinkling},\ref{subsec:both-wrinkling}), the sheet's size (through the parameter ${\cal R}$, Subsec. \ref{subsec:geo-limit}), and the boundary conditions at the far edge (Sec.~\ref{sec:Expansion}). Central to  all of these cases is the presence of a purely tensile, unwrinkled core, $0<r<L_I$, around the indenter, where the deformation is described by solving
the axisymmetric FvK equations (\ref{eq:FvK-nondim-2},\ref{eq:FvK-nondim-1}), subject to Eq.~(\ref{eq:wrinkled-I}) and the BCs~(\ref{eq:BC-sliding-wrinkling-I}), that yield two equations (\ref{eq:twomore}) for the three unknowns $\Psi(1),\tilde{a}, \tilde{L}_{I}$. The various cases in Subsecs.~\ref{subsec:II-wrinkling},\ref{subsec:both-wrinkling},\ref{subsec:geo-limit}, and Sec.~\ref{sec:Expansion} differ only in the final equation that connects $\Psi(1),\tilde{a}, \tilde{L}_{I}$, which stems from the continuity of radial displacement at the hole's edge (Eqs.~\ref{eq:thirdmatching}, \ref{eq:thirdmatching-1},\ref{eq:thirdmatching-1a},\ref{eq:thirdmatching-2}, respectively). In the following, we obtain the first two algebraic equations for $\Psi(1),\tilde{a}, \tilde{L}_{I}$, that are common to all of these cases.

Following Ref.~\cite{Vella15} (Sec.~3 of {\emph{Supplementary information}}), it is convenient to replace the dimensionless variables (\ref{eq:dimensionless-var}) with:
\begin{gather}
\bar{\rho} = \frac{r}{L_I} = \frac{\rho}{\tilde{L}_{I}}  \ \ ; \ \ \bar{\Psi} = \frac{\psi}{\sigma_{rr}(R)R} = \frac{\Psi}{\Psi(1)} \ \ ; \nonumber \\
\bar{\zeta} = \frac{z}{\sqrt{RL}} \sqrt{\frac{Y}{\sigma_{rr}(R)}}  =
\frac{\zeta}{\sqrt{\tilde{L}_{I} \Psi(1)}}
\label{eq:dimensionless-var-new}
\end{gather}
and the dimensionless force
${\cal F}$ (\ref{eq:dimensionless-param}) with:
\begin{equation}
\bar{\cal F}
=  \frac{1}{2 \pi R}
\sqrt{\frac{L_I}{R}} \sqrt{\frac{Y}{\sigma_{rr}(R)^3}} \ F   \ = \
{\cal F}\sqrt{\frac{\tilde{L}_{I}}{\Psi(1)^3}}  \ ,
 \label{eq:dimensionless-param-new}
\end{equation}
such that the BCs in (\ref{eq:BC-sliding-wrinkling-I},\ref{eq:twomore}) that involve explicitly the function $\bar{\Psi}(\bar{\rho})$ are:
\begin{equation}
\bar{\Psi}(\bar{\rho}=1) =1 \ \ ; \ \ \bar{\Psi}(\bar{\rho}=0)=0 \ \ ; \ \  \bar{\Psi}'(\bar{\rho}=1)=0 \ .
\label{eq:BC-new-new}
\end{equation}
Using a similar manipulation to the one employed earlier,
we make the additional transformation:
\begin{gather}
\Phi  = \bar{\rho} \bar{\Psi} \ \ ; \ \  \eta = \bar{\rho}^2 \ ,
\end{gather}
with which the BCs~(\ref{eq:BC-new-new}) become $\Phi(\eta = 0)=0 \ , \ \Phi(\eta = 1)=1 \ , \ \Phi'(\eta = 1)=\tfrac{1}{2}$, and the implicit expression for $\Phi(\eta)$, Eq.~(\ref{eq:FvK-nondim-tran-integrate-22}), is fully satisfied by the
numerical constants $A, \bar{\cal F}$, through the algebraic equations (\ref{eq:FvK-nondim-tran-integrate-2-edge},\ref{eq:FvK-nondim-tran-integrate-22-edge}). Solution of these equations yield the numerical values:
\begin{equation}
A \approx -0.697 \ \ ; \ \ \bar{\cal F} \approx 1.815 \ , \label{eq:nn11}
\end{equation}
which were found already in \cite{Vella15}. Equation~(\ref{eq:FvK-nondim-tran-integrate-1}), with $\zeta\to \bar{\zeta}$, together with the BCs for $\zeta$ in Eq.~(\ref{eq:BC-sliding-wrinkling-I}) yield:
\begin{gather}
\left(\tilde{a}\cdot (\tilde{L}_{I}-1) + \tilde{\delta}\right)\cdot \left(\tilde{L}_{I} \cdot \Psi(1)\right)^{-1/2} \nonumber \\
=\frac{2}{\sqrt{A}} \sinh^{-1}(\sqrt{A}) \approx 2.367 \ ,
\label{eq:new-new-11a}
\end{gather}
and the BC for the slope (\ref{eq:twomore}) becomes:
\begin{equation}
\tilde{a}\cdot \left(\tilde{L}_{I}/\Psi(1)\right)^{1/2} =  (1+A)^{-1/2} \approx 1.815  \ .
\label{eq:new-new-11b}
\end{equation}
For any value of the control parameter $\tilde{\delta} \gtrsim 3.3$, Eqs.~(\ref{eq:new-new-11a},\ref{eq:new-new-11b}), together with Eq.~(\ref{eq:thirdmatching}) for Subsec.~\ref{subsec:II-wrinkling}, or Eq.~(\ref{eq:thirdmatching-1}) for Subsec.~\ref{subsec:both-wrinkling}, or Eq.~(\ref{eq:thirdmatching-1a}) for Subsec.~\ref{subsec:geo-limit}, or Eq.~(\ref{eq:thirdmatching-2}) for Sec.~\ref{sec:Expansion}, form a set of 3 nonlinear algebraic equations for the three unknowns, $\Psi(1),\tilde{L}_{I},\tilde{a}$. The solutions of these equations
fully characterize the shape and stress of the deformed sheet in each case, and the corresponding indentation force is obtained with the aid of Eqs.~(\ref{eq:dimensionless-param-new},\ref{eq:nn11}).

\end{appendix}