\section{Why star clusters as cosmic ray sources? I -- Cosmic ray composition}

Differences between the composition of cosmic rays (CRs) and of the Solar System can be used to characterise the acceleration sites of energetic cosmic particles \cite{wiedenbeckreview,lingenfelterreview}.
Refractory elements are overabundant in CRs \cite{meyer} when compared to volatiles, requiring a preferential acceleration of grain material over diffuse gas (see the pioneering works \cite{epstein,cesarsky}, where two distinct acceleration mechanisms for refractories were described).
A common feature of both scenarios \cite{epstein,cesarsky} is that refractories are accelerated at supernova remnant (SNR) shocks, which are believed to be the main sources of the bulk of Galactic CRs (see \cite{hillas,me} for critical reviews of this paradigm).
The scenario proposed by Epstein \cite{epstein}, and later consolidated by Meyer, Drury, and Ellison \cite{meyer, ellison}, provides a more viable explanation of the overabundance of CR refractories, and will be briefly described here.

Dust grains are characterised by a very large mass-to-charge ratio, and therefore by a very large rigidity.
The correspondingly large Larmor radius makes dust grains easily injected and accelerated at SNR shocks \cite{epstein,meyer,ellison}.
Once accelerated, dust grains are eroded by sputtering.
Sputtered ions (i.e. refractory elements) are expected to move at the same velocity of the parent grain, which is likely to be suprathermal.
Therefore, sputtered refractories produced in the upstream region of a SNR will be advected back to the shock and, due to their large Larmor radius, further accelerated.
It follows that the amount of interstellar refractories which undergo acceleration and eventually become CRs is {\it not} (or not only) determined by the injection mechanism into the acceleration process (as it is quite plausible that most if not all of the grains swept up by the shock undergo acceleration), but rather by the effectiveness of sputtering.
Meyer and co-workers \cite{meyer,ellison} provided a very convincing argument that this scenario can explain the observed overabundance of refractories in CRs.
They even went further to claim that the relative abundance of CR volatile elements might be explained by an acceleration mechanism whose efficiency scales as the element mass-over-charge ratio.
Recent simulations showed that diffusive shock acceleration might indeed satisfy this requirement \cite{damiano}.
The outcome of Meyer, Drury and Ellison's work \cite{meyer} is that the composition of CRs is controlled by {\it both} volatility and mass-to-charge ratio.

A limitation of these pioneering approaches is that they were based on near-Earth observations of CRs (the only ones available at that time) which are heavily affected by Solar modulation at low (GeV domain and below) particle energies.
As a consequence, measurements of the composition of CRs were quite uncertain.
This problem was solved when the two Voyager probes crossed the heliopause, providing us with a measurement of the pristine spectrum of interstellar CRs (or, in case some residual modulation still persists, with something very close to it) down to particle energies of few MeV/nucleon~\cite{voyager,me2}\footnote{A puzzling outcome of CR transport models in the Galaxy is that, in order to explain CR data at both low (Voyager) and high (AMS-02 and other detectors) energies, CR sources must inject {\it broken power laws} in momentum, rather than the simple power laws $\propto p^{-s}$ predicted by shock acceleration theory \cite{voyager,vincentreview}.
It is convenient to express the CR injection spectrum for the element $i$ as a function of the particle energy per nucleon, $\dot{q}_i(E) \propto \beta^{-1} p^{-s}$.
For energies above the break $s \sim 4.2-4.3$ for all (primary) elements, while below the break $s \sim 4$ for $i =$~H and He, but it is much harder for heavies (e.g. $s \sim 3.3$ for $i =$~O).
Moreover, the break energy is in the multi-GeV domain for H and at few hundreds MeV/nucleon for heavier elements \cite{vincent}.
The origin of these different behaviours is not understood, and calls for further investigations.}.
These new accurate measurements of CR composition can be used to characterise the phase of the interstellar medium out of which CRs are accelerated \cite{eichmann,vincent}.

To understand how this can be done, let us follow Meyer and collaborators \cite{meyer,ellison} and suppose that an element of mass number $A$ and charge number $Z$ is accelerated with an efficiency scaling as the mass-to-charge ratio to some power: $\propto (A/Z)^{\alpha}$.
Simulations suggest, with quite some uncertainty, values of $\alpha = 1$ or 2 \cite{damiano}.
Consider then two extreme and idealised conditions characterising the gas reservoir out of which CRs are accelerated: an hot phase where all elements are fully ionised and a warm phase where all elements are singly ionised ($Z = 1$).
In the former case, all elements would be accelerated with roughly the same efficiency, as $A/Z \sim 2$, while in the latter the efficiency would scale with the mass number only, as $\propto A^{\alpha}$.
This illustrates how the ionisation state of the parent (pre acceleration) plasma impacts on the relative abundance of CR volatile elements.

The ionisation state is in turn set mostly by the (upstream) plasma temperature \cite{eichmann,vincent}, but a correction is needed if the acceleration operates at SNR shocks expanding in the warm phase of the interstellar medium.
This is because additional photoionisation of the ambient gas has to be taken into account.
It is due to UV and X-ray radiation from the supernova itself and from both the shocked ejecta and ambient medium (e.g. \cite{vincent} and references therein).
Once this correction is taken into account, the observed abundances of CR volatiles can be used to constrain the plasma temperature.
Tatischeff and collaborators \cite{vincent} found that data are best explained if the parent plasma is hot, with a temperature of the order of millions of degrees.
Such large temperatures are found inside superbubbles, i.e. the large cavities inflated in the interstellar medium due to the combined effect of stellar winds and recurrent supernova explosions in star clusters.
Moreover, the overabundance of CR refractories requires that such elements must be accelerated with an efficiency a factor of $\approx 10-50$ larger than that of volatiles.
To summarise, the results from Tatischeff and collaborators establish a clear link between volatile CRs and star clusters.
While the origin of refractories in CRs cannot be constrained in the absence of a clear knowledge of dust abundance and composition in the various phases of the interstellar medium, a scenario where {\it all} CRs (volatiles and refractories) are mainly accelerated in or around around star clusters is, to date, quite plausible (see discussion in \cite{vincent}).

While the scenario described above can reproduce the abundances of CR elements, it does not address the problem of {\it isotopic anomalies} which characterise energetic cosmic particles.
The most notable anomaly is the excess in the isotopic ratio $^{22}$Ne/$^{20}$Ne, which is equal to $\approx 0.3-0.4$ at CR sources, i.e. a factor of $\sim 4-5$ larger than in the Solar System, where the observed ratio is $\lesssim 0.1$\cite{binns,boschini}.
This prompted Casse and Paul \cite{casse} to suggest that a fraction of the CRs observed at the Earth must originate from chemically enriched wind material ejected by Wolf-Rayet (WR) stars.
As $^{22}$Ne is highly enhanced in WR stellar winds, some level of mixing of enriched (wind) and standard (Solar) material is required to reproduce the CR source isotopic ratio.

Three scenarios have been proposed to explain the $^{22}$Ne excess.
First, WR wind material can be accelerated by the SNR shock that forms after the star explodes \cite{meyer,prantzos}.
However, the model suffers from some fine tuning and, beside that, it is not clear whether WR stars eventually explode as supernovae or simply collapse into a black hole (see discussions and references in \cite{vincentreview,vincent}).
Second, the diffuse gas inside superbubbles may be chemically enriched due to the presence of WR stars in the parent cluster. Such material is then accelerated by SNR shocks as stars explode \cite{lingenfelterreview,higdon,binns}.
In this scenario, the composition of the material inside superbubble should be definitely non-solar, a thing which is not supported by X-ray observations of such objects (see e.g. \cite{jaskot}).
Finally, in the third and most promising scenario, enriched WR wind material is accelerated at the stellar wind termination shocks \cite{gupta,vincent}.
This latter possibility has been investigated quantitatively by Tatischeff and collaborators \cite{vincent}, who found that the $^{22}$Ne excess can be reproduced provided that $\approx 6$\% of the gas reservoir out of which CRs are accelerated is made of WR wind material.

Summarising, a link between volatile CRs and superbubbles (and therefore star clusters, which inflate such cavities) emerged from an accurate analysis of CR composition data.
Moreover, the most pronounced isotopic anomaly in the composition of CRs -- the excess in the $^{22}$Ne/$^{20}$Ne ratio -- can be explained by the acceleration of chemically enriched material at WR stellar wind termination shocks.
As massive stars are often found in clusters, it follows that also this isotopic anomaly is likely to originate in star clusters/superbubbles.

While it is not possible, at the moment, to decide whether also refractory CRs are originated at star clusters or not, the detection of 15 $^{60}$Fe nuclei in CRs might give us some hints \cite{60Fe}.
The $^{60}$Fe isotope is a short-lived radionuclide of half-life of $\sim$~2.6 Myr which is synthesised in core-collapse supernovae.
The fact that some energetic nuclei of this isotope reach the Earth means that their acceleration must have happened within the past few million years in the vicinity of the Solar system.
Assuming a {\it conventional} CR diffusion coefficient, it is possible to constrain the position of the acceleration site within a distance of $\approx$~1~kpc from the Sun (this is the displacement CRs would experience in a $^{60}$Fe half-life time).
A natural candidate for the acceleration of iron nuclei is the very nearby (less than 150 pc) Scorpius-Centaurus OB association \cite{60Fe}.
Recurrent supernova explosions took place in the association over the last few million years.
Such supernovae and their remnants might have been responsible for both the nucleosynthesis of $^{60}$Fe and its acceleration.
Moreover, the mechanical energy injected by stellar winds and supernovae in the Scorpious-Centaurus association is believed to be responsible for the formation of the Local Bubble within which the Solar system is located \cite{dieter}.
It is possible, even though far from being proven, that {\it all} CRs might come from star clusters/superbubbles, and we might even be located inside one of the accelerators of CRs!

\section{Why star clusters as cosmic ray sources? II -- Gamma-ray emission}

Enhanced gamma-ray emission is expected from CR sources, due to the decay of neutral pions produced when accelerated CRs interact with ambient matter.
In a pioneer paper, Montmerle~\cite{snob} used high energy gamma-ray observations to suggest that supernovae/OB associations (SNOBs) could be major sources of Galactic CRs.
This hypothesis remained for a long time much less popular than the mainstream scenario of CR acceleration at isolated SNRs, despite the fact that the number of associations of star clusters and superbubbles with GeV and TeV gamma-ray sources continued to grow slowly but steadily \cite{gammaclusters}.

The interest in stellar clusters as CR factories was revived when Aharonian and co-workers \cite{felixwinds}, after performing a morphological analysis of the high energy gamma-ray emission surrounding three of such objects, concluded that the emission had to be produced by CRs that had escaped the parent cluster after being accelerated there.
They also noticed that the gamma-ray spectrum of these objects extends to the multi-TeV domain without any sign of attenuation, implying that the spectrum  of the population of CRs responsible for the gamma-ray emission may reach the PeV energy domain and possibly beyond.
Finally, from the knowledge of the spatial distribution of the gas density, they could estimate the energy budget of CR in the regions, and concluded that the ensemble of {\it young}\footnote{The three objects they considered all have an age smaller than several million years.} star clusters in the Milky Way might provide a dominant contribution to the flux of the highest energy Galactic CRs.
As we will see in the next Section, young and old star cluster are expected to accelerate particles in a radically different way, and therefore it is important to increase the number of detections in gamma rays of objects of all ages.

In young clusters, the only source of mechanical energy are stellar winds.
It is a well known fact that the integrated power of all stellar winds in the Galaxy is a sizeable but most likely minor fraction of the power injected by supernova explosions \cite{winds}.
Therefore, a scenario could be envisaged where the bulk of Galactic CRs is produced at SNR shocks, while stellar winds provide the energy to accelerate the particles that we observe at the high energy end of the Galactic CR spectrum.
A way to test this scenario would be to observe in gamma-rays very young clusters (with an age smaller than very few Myr), to avoid any contamination from supernova explosions.
An attempt to do so can be found in \cite{giada}.

In very recent times, the number of star clusters and superbubbles detected in GeV gamma-rays has further increased \cite{gammaclustersnew}.
In the TeV domain, observations exploiting the superior angular resolution of imaging atmospheric Cherenkov telescope arrays revealed a complex and quite extended morphology around the star cluster Westerlund 1 \cite{westerlund1}.
Further high angular resolution studies will certainly help in identifying the sites of particle acceleration and constrain the acceleration mechanism.
Remarkably, 2 out of the 12 sources detected by the LHAASO detector at sub-PeV gamma-ray energies might be associated with young massive star clusters \cite{lhaaso}, reinforcing belief that such objects might operate as the most extreme Galactic CR accelerators.

\section{Why star clusters as cosmic ray sources? III -- Ideal environment for particle acceleration}

In order to illustrate why the environment surrounding star clusters is ideal for particle acceleration it is convenient to analyse separately the case of young and old clusters.
Following \cite{varenna}, we call young clusters those where no supernova has exploded yet, and for this reason their mechanical energy output is uniquely due to stellar winds.
A cluster becomes old few million years after its formation, when supernovae begin to explode.

For young clusters, it is necessary to further distinguish between compact and loose ones \cite{gupta,varenna,thibault}.
In a compact cluster, stars are so close to each other that individual stellar outflows merge to forme a collective wind which is decelerated at an expanding and (in zero-order approximation) spherical Wind Termination Shock (WTS).
In loose clusters each star develops its own wind and WTS.
The theory of particle acceleration at WTSs was first developed for individual stars \cite{wts}, and much later generalised to the collective WTS that forms around compact clusters \cite{gupta,giovanni}.
Interestingly, while the WTSs around individual stars are invariably strong, the collective one around compact clusters are characterised by a finite value of the Mach number, which increases monotonically as the WTS expands \cite{gupta,varenna}.
This might have important observational consequences, as the theory of diffuse shock acceleration predicts that stronger shocks accelerate particles with a harder energy spectrum.
A comparison between theoretical predictions and observations (e.g. gamma-ray spectra) is needed in order to constrain models and validate our knowledge on particle acceleration at WTS.
Note that the vast literature devoted to the theoretical description of particle acceleration at the Solar WTS might be of great help, due to the similarity of the problem \cite{modulation}.

Let us now discuss old clusters.
The environment surrounding such objects is a low density cavity (a superbubble) inside which recurrent supernova explosions take place.
As more and more stars explode, the mechanical energy input from stellar winds decreases accordingly.
The interactions between stellar winds and SNR shocks makes superbubbles highly supersonic and turbulent environments.
Theoretical descriptions of the acceleration of particles in superbubbles have been developed over the years, mainly by Andrii Bykov and co-workers \cite{andrii}, by other authors \cite{others} and, very recently, by Vieu and collaborators \cite{thibault,thibaulthillas,thibaultnew}.
The results of such theoretical investigations show that the expected spectra of accelerated particles are far from being universal and, depending on the richness of the cluster, acceleration might operate in an intermittent way.
This implies that a case by case study of individual clusters is certainly not sufficient, and that observations of many of them are needed in order to figure out how acceleration proceeds in superbubbles.

Another important issue is that of the maximum energy a particle can achieve when it undergoes acceleration around a star cluster.
A model independent estimate is provided by the Hillas criterium~\cite{hillas}, that states that the maximum energy $E_{max}$ scales with the size of the accelerator $R$, the strength of the magnetic field there $B$, and the characteristic fluid speed of the system $u$, or: $E_{max} \sim R u B$.
For a young (old) cluster, $R$ is the radius of the WTS (the size of the superbubble), and $u$ is the wind termination velocity (the velocity of turbulent motions inside the superbubble).
The magnetic field strength $B$ is often unknown, but can be constrained to be at most a sizeable fraction of all available mechanical/thermal energy in the system.
When that is done one sees that, under most optimistic conditions, compact young stellar clusters can accelerate protons barely up to the PeV domain, while turbulent motions inside superbubbles are unable to do so \cite{varenna,giovanni,thibaulthillas}.

These estimates of $E_{max}$ are disappointing, as the acceleration of protons well beyond PeV energies is required to explain the smooth observed transition between Galactic and extragalactic CRs.
A way to overcome this problem is to speculate that $\gg$~PeV particles are accelerated at very fast SNR shocks expanding within the collective wind launched by a compact cluster \cite{thibaultnew} (the magnetic field is expected to be large in such an environment \cite{badmaev}).
The required shock speed would be of the order of $u \approx 30000$~km/s \cite{thibaultnew}, but this would in turn correspond to an explosion energy of $E_{SN} \approx 10^{52} (M_{ej}/M_{\odot}) (u/30000~{\rm km/s})^2$~erg, i.e. an extremely energetic explosion.
Very powerful explosions are very rare, and this might be a difficult problem to overcome \cite{pierre}.

\section{Conclusions and open questions}

After having discussed CR composition, radiative signatures of ongoing acceleration at star clusters, and acceleration mechanisms in such objects, it is tempting to advance the hypothesis that most, if not all Galactic CRs might be produced inside superbubbles.
Proving or disproving this hypothesis will be certainly one of the most crucial issues to be addressed in coming years.
Some of the main issues to be solved are: {\it i)} the puzzling spectral features observed in the composition of low energy CRs that are, to date, very poorly understood, {\it ii)} the (long standing!) difficulty of models in explaining how particles are accelerated up to the PeV domain and beyond, and {\it iii)} the identification of the acceleration site of Galactic CRs of different energies (isolated or clustered SNR? WTS? Turbulent bubbles?).

\textit{Acknowledgements:} SG thanks the organisers of the Symposium and acknowledges support from Agence Nationale de la Recherche (grant ANR-21-CE31-0028).
He also thanks V. Tatischeff, T. Vieu, and G. Peron for many enlightening discussions.


\section*{Physical Experiments}

The photograph of our custom designed experimental apparatus is presented in Fig.~\ref{fig:exptSetup}A1-A2. Two main components -- a 3D-printed cylinder assembly and two vertically positioned linear translation stages -- make up the set. The hyperelastic materials of the bottom layer (VHB 4910, 3M) and the top layer (VHB 4950, 3M) have double-side stickiness. Material properties are listed in \textbf{Materials and Methods}. As shown in Fig.~\ref{fig:exptSetup}A2, the bottom layer is placed on top of the cylindrical casing. The contacting part is pressed by a 3d-printed gripper in in Fig.~\ref{fig:exptSetup}B1-B2.
The assembly is then placed between two vertical linear translation stages and the stages are moved to a specific height corresponding to a speicific prestretch. As the stages are moved upward, the piston -- pushed by the height adjuster -- also moves up but the outer edge of the bottom layer remains gripped onto the cylindrical casing. This leads to a uniform biaxial stretch in the hyperelastic material. The relation between the imposed stretch and height can be determined from simple geometry. The top layer is glued onto the stretched bottom layer to form a bilayer composite. The bilayer composite is then cut out from the apparatus and released from the gripper. This resulting 3D structure is referred to as the ``soft kirigami composite" in this paper.

\section*{Soft Circular Composites}\label{sec:SoftCircularComposites}

Our investigation starts with the setup described in Fig.~\ref{fig:overview}. A bottom layer of circular shape is first stretched by a prescribed amount of prestretch ($\lambda=1$ corresponds to stress-free configuration). A stress-free circular top layer is affixed onto the bottom layer. Physical parameters are described in \textbf{Materials and Methods.} Once the composite structure is released from the experimental setup of Fig.~\ref{fig:exptSetup}, a variety of shapes emerge depending on prestretch, $\lambda$. In Fig.~\ref{fig:circularComposites}A, a series of four shapes of the experimental composites are presented. These shapes -- that only differ by $\lambda$ -- are qualitatively different. The number of wrinkles on the outer edge of the composites, represented by the $k$ number in the figure, increases from $k=0$ in Fig.~\ref{fig:circularComposites}A1 to $k=4$ in Fig.~\ref{fig:circularComposites}A4. Fig.~\ref{fig:circularComposites}B shows finite element (FE) simulations that can also capture the qualitatively distinct shapes. While developing FE simulations, we realized the existence of different branches in the solution space and a propensity for the FE method to reach a local energy minimum instead of the global. This exposes a challenge in numerical simulation of such structures and motivates us to fundamentally understand the problem through a theoretical lens. The next section outlines the theory that was developed to explain the observations in Fig.~\ref{fig:circularComposites}.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.98\textwidth]{figure3_ver12_Wenzhong.pdf}
    \caption{
    %
    Free buckling shapes of soft composites. (A) Experimental images and (B) snapshots from finite element simulations at different values of prestretch: (1) $\lambda$ = 1.015, (2) $\lambda$ = 1.064, (3) $\lambda$ = 1.108, and (4) $\lambda$ = 1.167 to obtain different number of waves ($k=0, 2, 3$, and $4$, respectively) on the outer edge. Scalebar: 1 cm.
    %
    For experimental parameters see Section~\ref{sec:parameters}.
    }
    \label{fig:circularComposites}
\end{figure}

\begin{figure}[h!]
\centering
   \includegraphics[width=0.90\textwidth]{Figures/Fig_Theory.pdf}
   \caption{
   Theory with branches.
   The bottom figure plots the $k$ value as a function of prestretch, $\lambda$. For a given geometric and material properties there exists a critical $k$ after which the only possible (minimal energy) solution corresponds to the scroll-mode (i.e., there exists a max $k$).
   }
   \label{fig:fig_theory}
\end{figure}

\section*{Composite plate theory}

We assume plane stress and Kirchhoff-Love (KL) kinematic assumption on total in-plane displacement $\doubleunderline{u}(r,\theta,z) = \doubleunderline{u}^0(r,\theta) - z \nabla w(r,\theta)$ as the circular laminated composite plates used in experiments are thin. Here, $r$, $\theta$, $z$ are coordinates of the plate in radial, circular and direction across the thickness, $\doubleunderline{u}^0$ and $w$ are the displacement vector of the mid-surface and its vertical displacement, respectively.
Furthermore, we use the Green-Lagrange strain tensor under F\"{o}ppl-von Karman (FvK) kinematic assumptions, which further simplifies the theory and limits it to small strains and moderate rotations,
\begin{equation}
\doubleunderline{E}_\text{FvK} = \doubleunderline{E}_\text{FvK}^0 - z \nabla \otimes \nabla w,
\end{equation}
where $2 \doubleunderline{E}_\text{FvK}^0 = \nabla \doubleunderline{u}^0 + \left ( \nabla \doubleunderline{u}^0 \right )^\mathrm{T} + \nabla w \otimes \nabla w$ is the GL strain tensor of the mid-surface under FvK assumptions. In these expressions, $\nabla \otimes \nabla w$ is the curvature tensor, which will be later denoted by $\doubleunderline{\kappa}$ and $\nabla w \otimes \nabla w$
is the only nonlinear term remaining from the complete GL strain tensor.
For the material we assume that it is compliant with St. Venant-Kirchhoff elastic material model,
\begin{equation}\label{eqn:constitutive_equation}
	\doubleunderline{\sigma}_{K} = \frac{E_K}{1-\nu_K^2}\left(
	(1-\nu_K) \left(\doubleunderline{E}_\text{FvK} - \doubleunderline{S}_{K} \right) + \nu_K \tr \left(\doubleunderline{E}_\text{FvK} - \doubleunderline{S}_{K} \right) \doubleunderline{I}
	\right),
\end{equation}
where $E_K$, $\nu_K$, $\doubleunderline{S}_{K}$ are Young modulus, Poisson ratio and prestrain tensor in $K$-th layer in the composite. As in experiment, we apply $\doubleunderline{S}_{K} = - \doubleunderline{I} (\lambda - 1)$, where $\lambda$ and $\doubleunderline{I}$ are the applied prestretch and a unity matrix, respectively.

A membrane force and bending moment (per unit length) tensors are as usually defined by $\doubleunderline{N} = \int_h \doubleunderline{\sigma} dz$ and $ \doubleunderline{M} = \int_h \doubleunderline{\sigma} z dz$, respectively, where $h$ is the overall thickness of the composite plate.
A coupling between between $\doubleunderline{N}$, $\doubleunderline{M}$ as functions of $\doubleunderline{\kappa}$, $\doubleunderline{E}_{\text{FvK}}^0$ and the applied prestretch $\lambda$ is obtained
\begin{equation}\label{eqn:NMsistem}
    \begin{aligned}
        \doubleunderline{N} = &\; \doubleunderline{\overline{\alpha}}
         \doubleunderline{E}_{\text{FvK}}^0 +
        \doubleunderline{\overline{\beta}} \doubleunderline{\kappa} + \doubleunderline{\overline{\gamma}}(\lambda),\\
        \doubleunderline{M} = & \doubleunderline{\tilde{\alpha}} \doubleunderline{E}_{\text{FvK}}^0 +
        \doubleunderline{\tilde{\beta}} \doubleunderline{\kappa} + \doubleunderline{\tilde{\gamma}}(\lambda).
    \end{aligned}
\end{equation}
Furthermore, expressing $\doubleunderline{E}_{\text{FvK}}^0$ from \eqref{eqn:NMsistem}$_1$ and inserting it into \eqref{eqn:NMsistem}$_2$ yields expression for $\doubleunderline{M}$ as a function of $\doubleunderline{\kappa}$ and $\doubleunderline{N}$,
\begin{equation}\label{eqn:Main_compact_form_short}
    \begin{aligned}
        \doubleunderline{M} = \doubleunderline{{e}} \doubleunderline{N} -
        \doubleunderline{\overline{D}} \doubleunderline{\kappa} + \doubleunderline{{\gamma}}_0.
    \end{aligned}
\end{equation}
Here, symbols $\doubleunderline{\overline{\alpha}}$, $\doubleunderline{\overline{\beta}}$, $\doubleunderline{\overline{\gamma}}$,
$\doubleunderline{\tilde{\alpha}}$, $\doubleunderline{\tilde{\beta}}$, $\doubleunderline{\tilde{\gamma}}$, $\doubleunderline{e}$, $\doubleunderline{\overline{D}}$, $\doubleunderline{\gamma}_0$ denote auxiliary functions that are defined later in this section, see Eqs. (\ref{eqn:stiffnessesTensors}) and (\ref{eqn:stiffnesses}).

Let us define a stress function $\mathcal{F}$ such that
\begin{equation}\label{eqn:nablaFN}
    \tilde{\nabla} \tilde{\nabla} \mathcal{F}  := \Delta \mathcal{F}\doubleunderline{I} - \nabla\otimes\nabla \mathcal{F} = \doubleunderline{N}
\end{equation}
and that satisfies the in-plane equilibrium
$\nabla \cdot \doubleunderline{N} = \doubleunderline{0}$.
Inserting Eqs. (\ref{eqn:Main_compact_form_short}) and (\ref{eqn:nablaFN}) into the out-of-plane equilibrium of the plate $\nabla \cdot (\nabla \cdot \doubleunderline{M}) + \nabla \cdot \left ( \doubleunderline{N} \cdot \nabla w\right) = 0$ its alternative form is obtained
\begin{equation}\label{eqn:Main_eqnsF}
    \overline{D} \Delta^2 w  - \nabla \cdot \left (\tilde{\nabla} \tilde{\nabla} \mathcal{F}  \nabla w\right)  =  0.
\end{equation}
Similarly, an alternative form of the compatibility conditions $\nabla \times \left( \nabla \times \doubleunderline{E}^0_{\text{FvK}} \right)^{\mathrm{T}} = \boldsymbol{0}$ can be written as
\begin{equation}\label{eqn:komp}
    \Delta^2 \mathcal{F} + \frac{\overline{\alpha} (1-\overline{\nu}_A^2)}{2} \left[ w , w\right] = 0,
\end{equation}
where $\left[\phantom{.},\phantom{.}\right]$ is the Monge-Ampere operator.

\begin{equation}\label{eqn:stiffnesses}
    \begin{aligned}
        \tilde{\alpha} & =  \left (  \frac{E_2}{1-\nu_2^2} e_2 h_2 - \frac{E_1}{1-\nu_1^2} e_1 h_1 \right )\\
        \tilde{\beta} & =  - \left (  \frac{E_2}{1-\nu_2^2} \left (\frac{h_2^3}{12}  + h_2 e_2^2\right ) + \frac{E_1}{1-\nu_1^2} \left (\frac{h_1^3}{12}   + h_1 e_1^2 \right ) \right )\\
        \tilde{\gamma} & =  - \left (  \frac{E_2}{1-\nu_2}  h_2 e_2 \eta_2 - \frac{E_1}{1-\nu_1}  h_1 e_1 \eta_1       \right )\\
        \tilde{\nu}_A & =  \frac{1}{\tilde{\alpha}}
                \left (  \frac{E_2 \nu_2}{1-\nu_2^2} e_2 h_2 - \frac{E_1 \nu_1}{1-\nu_1^2} e_1 h_1 \right )\\
        \tilde{\nu}_B & =  \frac{-1}{\tilde{\beta}} \left (  \frac{E_2 \nu_2}{1-\nu_2^2} \left (\frac{h_2^3}{12}  + h_2 e_2^2\right ) + \frac{E_1 \nu_1}{1-\nu_1^2} \left (\frac{h_1^3}{12}   + h_1 e_1^2 \right ) \right )\\
        \overline{\alpha} & =   \left (  \frac{E_2}{1-\nu_2^2} h_2 + \frac{E_1}{1-\nu_1^2}  h_1 \right )\\
        \overline{\beta} & =   - \tilde{\alpha}\\
        \overline{\gamma} & =  - \left (  \frac{E_2}{1-\nu_2}  h_2  \eta_2 + \frac{E_1}{1-\nu_1}  h_1  \eta_1       \right )\\
        \overline{\nu}_A & =  \frac{1}{\overline{\alpha}}
                \left (  \frac{E_2 \nu_2}{1-\nu_2^2}  h_2 + \frac{E_1 \nu_1}{1-\nu_1^2}  h_1 \right )\\
        \overline{\nu}_B & =  \tilde{\nu}_A
    \end{aligned}
\end{equation}

\begin{equation}\label{eqn:stiffnessesTensors}
    \begin{aligned}
        \doubleunderline{\overline{\alpha}} & =
        \overline{\alpha}
        \begin{bmatrix}
        1  & \overline{\nu}_A & 0\\
         \overline{\nu}_A  & 1  & 0 \\
        0 & 0 &  2(1 - \overline{\nu}_A)
        \end{bmatrix}\\
        \doubleunderline{\overline{\beta}} & =
        \overline{\beta}
        \begin{bmatrix}
        1  & \overline{\nu}_B & 0\\
         \overline{\nu}_B  & 1  & 0 \\
        0 & 0 &  2(1 - \overline{\nu}_B)
        \end{bmatrix}\\
        \overline{\doubleunderline{\gamma}} & = \doubleunderline{I}\; \overline{\gamma}\\
        \doubleunderline{\tilde{\alpha}} & =
        \tilde{\alpha}
        \begin{bmatrix}
        1  & \tilde{\nu}_A & 0\\
         \tilde{\nu}_A  & 1  & 0 \\
        0 & 0 &  2(1 - \tilde{\nu}_A)
        \end{bmatrix}\\
        \doubleunderline{\tilde{\beta}} & =
        \tilde{\beta}
        \begin{bmatrix}
        1  & \tilde{\nu}_B & 0\\
         \tilde{\nu}_B  & 1  & 0 \\
        0 & 0 &  2(1 - \tilde{\nu}_B)
        \end{bmatrix}\\
        \tilde{\doubleunderline{\gamma}} & = \doubleunderline{I}\; \tilde{\gamma}.
    \end{aligned}
\end{equation}

\subsection{Uncut symmetric solution}\label{subsec:sim_res}
Following experimental results shown in Fig. \ref{fig:circularComposites} we first seek an axi-symmetric solution, which we refer to as a ``cup''. The theory suggests that after the prestretched layer is adhered to a nonprestretched one the membrane stresses are equlibrated in both layers after the release, but due to the coupling between membrane stresses in the mid-surface of each layer, bending moments are induced.
These lead to the shortening and compressive (circular) membrane forces in the circumference of the two-layered circular composite plate, most prominently towards the edge. By plane equilibrium, tensile membrane forces are induced towards the center of the circular plate.

Due to the symmetry of the problem, the order of both PDEs (\ref{eqn:Main_eqnsF}) and (\ref{eqn:komp})  that describe the deformation of the composite plate can be reduced by substitution $\varphi = w_{,r}$ and $\phi = \mathcal{F}_{,r}$. After integration we obtain
\begin{equation}\label{eqn:Main_2nd_order_eq}
\begin{aligned}
    \overline{D}  r \left( \Delta \varphi - \frac{\varphi}{r^2} \right)
    -\phi \varphi &= 0 \\
    \frac{1}{\overline{\alpha}(1-\overline{\nu}_A^2)} r \left ( \Delta \phi - \frac{\phi}{r^2} \right) + \frac{\varphi^2}{2} &= 0.
 \end{aligned}
\end{equation}
We approximate $\varphi$ with the model function $\varphi(r) = C_1 r^n$, to obtain $\phi(r) = C_1^2  r/a ( 1 - ( r/a) ^{2 n}) / (8 \Lambda  n (n+1))$ from \eqref{eqn:Main_2nd_order_eq}$_2$.
From the boundary conditions $M_{rr} = 0$ and $N_{rr} = 0$ at the edge, we obtain $C_1 =  \gamma_0 a / (\overline{D}(\nu + n))$. Here, we took $\nu = \tilde{\nu}_A \doteq \tilde{\nu}_B \doteq \overline{\nu}_A \doteq \overline{\nu}_B$ for the sake of simplicity. The exponent $n$ can be determined by applying orthogonality condition on \eqref{eqn:Main_2nd_order_eq}$_1$,
$\int_{0}^{a} \left ( \overline{D} r \left ( \Delta \varphi - {\varphi}/{r^2} \right ) - 	\phi    \varphi \right ) {\partial \varphi}/{\partial n}   dr =  0$.
This condition yields an implicit relation between $n$, the dimensionless parameter $\chi$ and $\nu$ in the following form $3 (2 n+1)^2 \chi -4 \left(n^2-1\right) \left(8 n^2+18 n+9\right)^2 (\nu +n)^2 = 0$, where $\chi = B^2 a/(\Lambda \overline{D})$, $B = \gamma_0 a / \overline{D}$ and $\Lambda = 1/((1-\overline{\nu}_A^2)\overline{\alpha} a)$, while $\gamma_0$ is the moment due to prestretch. The explicit expression for $n(\chi(\lambda))$ is solved numerically.

From experiments and simulations described earlier we learn that by increasing the prestrecth, the axisymmetric cup structure losses its stability and wrinkles into a $k$-fold axisymmetric structure when the circular forces exceed some critical.

\subsection{\emph{k}-fold axisymmetric postbuckling}
We seek the $k$-fold axisymmetric solution with the following two model functions $w(r,\theta) = \overline{w}(r) + \tilde{w}(r,\theta)$ and $\mathcal{F}(r,\theta)  = \overline{\mathcal{F}}(r) + \tilde{\mathcal{F}}(r,\theta)$, where
$\overline{w}(r) = \int \varphi(r) dr$ and $\overline{\mathcal{F}}(r) = \int \phi(r) dr$ are the axisymmetric solutions from before and $\tilde{w}(r,\theta) = f(r) \cos k \theta$ and $\tilde{\mathcal{F}}(r,\theta) =  g(r) \cos k \theta$ are the $k$-fold axisymmetric eins\"{a}tze. Plugging both model functions into Eqs. (\ref{eqn:Main_eqnsF}) and (\ref{eqn:komp}) yields
\begin{equation}\label{eq:Main_MDMVTheory}
	\begin{aligned}
		\overline{D} \Delta^2 \tilde{w} -  \doubleunderline{\kappa}_0 : (\nabla \otimes \nabla{\tilde{\mathcal{F}}})
		- \doubleunderline{N}_0 : (\nabla \otimes \nabla \tilde{w} )
		-  \left[ \tilde{\mathcal{F}} , \tilde{w}  \right] & = 0, \\
		\frac{\Delta^2 \tilde{\mathcal{F}} }{\overline{\alpha}(1-\overline{\nu}_A^2)}  +
		\doubleunderline{\kappa}_0 : (\nabla \otimes \nabla \tilde{w})  +
		\frac{1}{2} \left[ \tilde{w}  , \tilde{w}  \right]
		& = 0.
	\end{aligned}
\end{equation}
Here, $\doubleunderline{\kappa}_0 = \Delta \overline{w} \doubleunderline{I} - \nabla \otimes \nabla \overline{w}$ and $\doubleunderline{N}_0$ are modified curvature tensor and membrane force tensor of the symmetric solution, respectively.

We observe that for smaller prestretches $\lambda$ the first most energetically favourable (stable) wrinkling mode is 2-fold symmetric (Pringles chip-like), followed by 3-fold, 4-fold etc.  symmetric solutions, when the prestretch is increased.
Fig. \ref{fig:fig_theory}, where the occuring deformation modes and their strain energies are displayed. With each following mode the wrinkles are more and more localized at the edge, because with deeper and deeper spherical cap, the membrane compressive stresses are more and more localized at the edge (see Fig. \ref{fig:InnerMembraneForces}) and then released through wrinkling. Because wrinkles with large $k$-fold symmetry are more localized at the edge, they are energetically more favorable, because they release more of the membrane energy, which is localized at the edge.

Another reason why the higher deformation modes $k$ are promoted with larger prestretches $\lambda$ is that the increasing curvature of the system (spherical cup) provides geometric rigidity, which acts analogously to an elastic substrate. This is evident if we linearize \eqref{eq:Main_MDMVTheory}$_2$ and insert it into \eqref{eq:Main_MDMVTheory}$_1$. We obtain an equation of the form $\tilde{D} \Delta^2 \tilde{w} - \doubleunderline{N_0} : \nabla \otimes \nabla \tilde{w} + \alpha_{0} (\kappa_0:\kappa_0) \tilde{w} = 0$, which is analogous to the equation that governs wrinkling of thin compressed films on compliant substrates. The therm $\alpha_{0} (\kappa_0:\kappa_0) \tilde{w}$ penalizes buckling modes with large amplitude to wavelength ratios, therefore buckling modes with larger $k$ are energetically preferable.

\subsection{Deep postbuckling}
When the stretching $\lambda$ of the membrane is large we observe that $k$-fold symmetric solution is no longer energetically favorable. Instead, we detect experimentally and computationally (using Abaqus) that cylindrical bending solution is dominant. We refer to this solution as a ``scroll'' solution. This reduces Eq.~\eqref{eqn:Main_eqnsF} to
$\kappa_{11,11}  =  0$, as only $\kappa_{11} \neq 0$ and $\kappa_{22} = \kappa_{12} = 0$, while Eq.~\eqref{eqn:komp} is automatically satisfied.
Assuming stress free edges, we obtain $\doubleunderline{N} = \doubleunderline{0}$ and $\kappa_{11} = \gamma_0/\overline{D}$. The strain energy is calculated as $\mathcal{E} = \pi a^2 \gamma_0^2/(2\overline{D})$.

The $k$-fold symmetries do not increase indefinitely. At a certain point the scroll solution (which is similar to the 2-fold symmetry) offers the largest relaxation of the membrane stresses. Its membrane stresses are released entirely, however the price for this is a bit larger bending strain energy, than in the $k$-fold symmetric cases, therefore it is not preferable with smaller $\lambda$.

\begin{figure*}[h!]
\centering
    \includegraphics[width=1\textwidth]{Figure5_Kirigami_arxiv.pdf}
    \caption{
    Soft kirigami composites. (A) Experimental images and (B) snapshots from finite element simulations at different values of prestretch: (1) $\lambda = 1.05$, (2) $\lambda=1.15$, (3) $\lambda=1.2$, and (4) $\lambda=1.25$ to obtain the preprogrammed shapes. It turns out that a (quasi)axisymmetrical shape is not energetically favorable for $1.15 \lesssim \lambda \lesssim 1.25$ for this particular set of material and geometric parameters.
    (C) Normalized height, $\bar H = H/R_o$, as a function of prestretch, $\lambda$, where $R_o$ is the radius of the kirigami layer.
    %
    (D) Normalized strain energy, $\bar E$, as a function of $\lambda$ at two different mode shapes: $k=0$ and $k=2$.
    The strain energy has been normalized by the characteristic bending energy, as described in the text. The inset shows the difference in normalized strain energy between the two ($\Delta \bar E = \bar E_{k=0} - \bar E_{k=2}$, where the subscript indicates the mode number) as a function of $\lambda$.
    %
    (E) Critical prestretch, $\lambda^{k=2}_{\textrm{crit}}$, as a function of the portion of the removed material, $\alpha$ ($\alpha=1$ means the entire layer has been removed).
    }
    \label{fig:kirigamiComposites}
\end{figure*}

\section*{Soft Kirigami Composites}

Informed by the nonlinear nature of the problem and the existence of multiple branches in the previous section, we essentially need to find out a way of reaching the hemispherical cap solution. Therefore, we need to achieve at least ``quasi-axisymmetric", if not an axi-symmetric ($k=0$) solution at at larger values of prestretch. A opposed to the 1D example, where applying enough prestretch enables us to produce a circle, here on the circular plate circular compressive forces \cite{stein2019buckling, Holmes2019COCIS, Pezzulla2015SM} occuring due to the change of circumference during bending, prevent that the circular plate would bend into a hemisphere, as the $k$-fold wrinkling inastability or scroll instability occurs at the edge \cite{Hadrien2017SM, Holmes2019COCIS,Pezzulla2015SM}.

The first solution that comes to mind \cite{moshe2019kirigami, tang2019programmable, cui2018origami} is therefore to cut the top layer along a prescribed path, whereby forming a so called ``soft kirigami composite" as a way of releasing the compressive force from certain areas of the composite and partially remove the wrinkled nature of the outer edge.
A heuristic shape for the top layer was chosen as shown in Fig.~\ref{fig:overview}B.
Physical parameters are presented in \textbf{Materials and Methods}. Fig.~\ref{fig:kirigamiComposites}A and B present snapshots of the structural shapes at different values of prestretch from experiments and simulations, respectively.

If we remove a circumferential fraction of material $\alpha$, which is divided between $n_c$ pieces, we enforce $n_c$-fold symmetry, which is unfortunatelly only  ``quasi-axisymmetric" and not axi-symmetric $k=0$.  This is because, even at small prestretches $\lambda$, but even more so at the prestretches larger than the critical prestretch for $n_c$-fold buckling $\lambda > \lambda_\text{crit}^{k=n_c}$ the inhomoegnity in the circular direction causes additional bending of the ``pizza slices''  in the circular direction. Therefore to obtain a solution that is as axi-symmetric as possible a large number $n_c \to \infty$ pieces of composite have to be removed.
However this might not be the best idea since both, the number of cuts $n_c$ and the choice of the cut out ratio $\alpha$ affect weather the desired $n_c$-fold symmetric hemisphere is stable or weather it looses stability and buckles into an undesired low $k$-fold symmetric deformation mode.
Therefore we are searching for the appropriate $\alpha$ and $n_c$  to ensure that:
 1.) That the $n_c$ -fold symmetric hemisphere is stable.
 2.) That it is symmetric enough, since 2,3,4,5-fold symmetric shapes are not symmetric enough for our liking.

We employed an empirical trial and error approach to find a suitable set of parameters for a shape close to a hemispherical cap in Fig.~\ref{fig:kirigamiComposites}A-B at $\lambda =1.25$. Fig.~\ref{fig:kirigamiComposites}C and D show the height, $H$, of the kirigami composite and the elastic energy, $E$, as functions of prestretch, $\lambda$. The height has been normalized by the radius of the substrate ($R=2$ cm) so that $\bar H = H / R$. The strain energy has been normalized by the characteristic elastic energy of the composite as discussed in \textbf{Materials and Methods.} The difference in energy between two different modes ($k=0$ and $k=2$) is too small to be visually observable, therefore, Fig.~\ref{fig:kirigamiComposites} includes an inset showing the normalized difference in strain energy between the two modes: $\Delta \bar E = \bar E_{k=0} - \bar E_{k=2}$. Note that $\Delta \bar E$ is positive between $1.15 \lesssim \lambda \lesssim 1.30$; i.e., $k=2$ is energetically favorable compared with $k=0$ shape.

This critical value of prestretch, $\lambda_\textrm{crit}^{k=2}$, at which $k=2$ shape appears from a $k=0$ shape depends on the amount of removed material in the Kirigami layer, $\alpha$, among other parameters. $\alpha=1$ means the entire layer has been removed, whereas $\alpha=0$ corresponds to a circular shape with no Kirigami cuts at all. Fig.~\ref{fig:kirigamiComposites}E shows the critical prestretch, $\lambda_\textrm{crit}^{k=2}$, as a function of $\alpha$ from finite element simulations and theoretical analysis.

\section{Materials and Methods}

\subsection*{Physical Parameters}
A commercially available acrylic adhesive (3M VHB tape) is used as the bottom and kirigami layers. Mooney-Rivlin model is used in the finite element simulations of Fig.~\ref{fig:circularComposites} and Fig.~\ref{fig:kirigamiComposites}. The strain energy per unit of reference volume in this model is
\begin{equation}
    U = C_{1} \left( \bar I_1 - 3 \right) + C_{2} \left( \bar I_2 - 3 \right) + D_1 \left( J - 1 \right)^2,
\end{equation}
where $C_1, C_2$, and $D_1$ are material parameters, $\bar I_1$ and $\bar I_2$ are the first and second deviatoric strain invariants, and $J$ is the elastic volume ratio. The bottom layer has a thickness of $h=1.0$~mm and its material parameters are $C_1 =4.84$~kPa, $C_2=14.536$~kPa, and $D_1 = 0.96\times10^{3}$~kPa. The kirigami/top layer is $h=1.1$~mm thick with the following material parameters:
$C_1 = 4.00$~kPa, $C_2 = 78.7152$~kPa, and $D_1=4.1082\times10^{3}$~kPa.

In order to normalize the elastic strain energy, we formulate a  characteristic elastic energy of the composites: $E^* = Y h^3$, where $Y = 2 G (1 + \nu)$ is the effective Young's modulus of the substrate, $h$ is the thickness of the substrate, and $\nu=\frac{3K-2G}{6K+2G}$, $G = 2(C_1 + C_2)$, and $K=2/D1.$

For the soft kirigami composites in Fig.~\ref{fig:kirigamiComposites}, the number of cuts is  $n_c=10$ and the ratio of the inner radius to the outer radius is $\bar r = 0.4$ and $d\theta=0.1$ radians, i.e., $1$ out of $2\pi$ radian was removed. This corresponds to $\alpha=0.13$.