# Non-degenerate surface pair density wave in the Kagome superconductor
CsV3Sb5 \- application to vestigial orders

Yue Yu  Department of Physics, Stanford University, Stanford, CA 94305
Department of Physics, University of Wisconsin, Milwaukee, WI 53201

###### Abstract

On the Sb-layer of the Kagome superconductor CsV3Sb5, pair density wave states
have been observed. When the high-temperature charge orderings are treated as
static backgrounds, these PDW states exhibit the same wavevector in the
effective 2D Brillouin zone. Interestingly, these PDW states break the same
symmetry on the surface. Considering the presence of this non-degenerate PDW,
we investigate the implications for the possible existence of a vestigial
charge-4e phase with a non-zero center-of-mass momentum. To distinguish
between different vestigial phases, we propose scanning tunneling microscopy
experiments. We aim to provide insights into the nature of the vestigial
phases and their distinct characteristics in CsV3Sb5. This research sheds
light on the interplay between PDW states, charge orderings, and
superconductivity of the Kagome superconductor.

##  I  introduction

The recently discovered Kagome superconductor CsV3Sb5 (CVS) Ortiz et al.
(2019, 2020) is a highly intriguing material that has attracted significant
experimental and theoretical interest due to its exotic charge and
superconducting orderings Jiang et al. (2023).

In its high-temperature phase, the CVS crystalizes in the P6/mmm space group
and exhibits a layered structure composed of V-Sb sheets intercalated by Cs
atoms Ortiz et al. (2019). At a temperature of 94K, scanning tunneling
microscopy (STM) experiments have revealed a $2a_{0}$ charge density wave
(CDW) state with a $2\times 2$ superlattice modulation Xu et al. (2021); Chen
et al. (2021a); Liang et al. (2021); Wang et al. (2021a). Notably, the CDW has
been found to possess a three-dimensional character, with reports of both
$2\times 2\times 2$ and $2\times 2\times 4$ modulations Li et al. (2021); Miao
et al. (2021); Ortiz et al. (2021a); Liang et al. (2021); Song et al. (2022).

Moreover, different intensity distributions of the CDW peaks with unusual
magnetic response have been observed, indicating the presence of a chiral
charge order Jiang et al. (2021); Zhao et al. (2021a); Wang et al. (2021a).
Investigations are underway on the time-reversal symmetry breaking in this
system Mielke III et al. (2022); Yu et al. (2021a); Li et al. (2022a); Feng et
al. (2021a); Park et al. (2021); Feng et al. (2021b); Denner et al. (2021);
Lin and Nandkishore (2021); Yu et al. (2021b). The difference in CDW intensity
also suggests rotational symmetry breaking. The temperature dependence of the
nematic transition has been recently explored Xiang et al. (2021); Ni et al.
(2021); Li et al. (2021); Ratcliff et al. (2021); Wang et al. (2021b); Luo et
al. (2022a).

Around 60K, an unidirectional CDW with a periodicity of $4a_{0}$ is observed
on the Sb-layer through STM measurements Zhao et al. (2021a); Liang et al.
(2021); Chen et al. (2021a); Li et al. (2023). However, this CDW is not
observed on the alkali layer or in the bulk Li et al. (2021); Miao et al.
(2021); Ortiz et al. (2021a); Song et al. (2022). Superconductivity is
observed in CsV3Sb5 at a temperature of 2.3K Ortiz et al. (2020), and the
nature of its order parameter is still under investigationMu et al. (2021); Yu
and Li (2012); Kiesel et al. (2013); Wang et al. (2013); Duan et al. (2021);
Luo et al. (2022b); Xu et al. (2021); Zhao et al. (2021b); Wu et al. (2021);
Chen et al. (2021a); Gu et al. (2022); Gupta et al. (2022); Chen et al.
(2021b); Yu et al. (2021c); Chen et al. (2021c); Zhang et al. (2021).

In the context of the superconducting phase, the presence of pair density wave
(PDW) on the Sb-layer has been reported by STM probes Chen et al. (2021a).
Similar to the $4a_{0}$ CDW, its existence on the alkali layer and in the bulk
remains unconfirmed, resulting in unknown 3D wave vectors. Consequently, a
comprehensive analysis of the bulk superconducting phase diagram becomes
challenging. Conversely, the STM results offer valuable insights and provide
sufficient information to construct the superconducting phase diagram on the
surface. Therefore, studying the surface phase diagram can serve as a stepping
stone towards a better understanding and analysis of the bulk superconducting
phase.

In this study, our focus is directed towards the surface superconducting
transitions. To simplify the analysis, we consider the charge orderings that
emerges at significantly higher temperatures than the superconductivity as a
static background. This treatment effectively enlarges the unit cell and
results in a folded Brillouin zone. Within this folded Brillouin zone, it is
observed that all PDWs exhibit the same wavevector. These PDWs, being
physically equivalent, can be effectively described by a single order
parameter denoted as $\Delta_{\bf Q}$.

The presence of the unusual PDW has significant implications for the vestigial
phases. Traditionally, vestigial phases are constructed using two independent
PDW order parametersBerg et al. (2009). However, in the case of the non-
degenerate PDW $\Delta_{\bf Q}$ and the uniform SC $\Delta_{\bf 0}$ in CVS,
both order parameters need to be considered to construct vestigial phases.

Specifically, if a charge-4e phase exists, it would be characterized by the
composite order parameter $\Delta^{4e}=\Delta_{\bf 0}\Delta_{\bf Q}$. The
presence of such a phase can be examined through experimental techniques other
than the Little-Parks oscillationGe et al. (2022). In this study, we will
utilize STM signatures to differentiate between different vestigial phases and
provide insights into their distinct characteristics.

This paper will primarily investigate the properties of the non-degenerate PDW
state on the low-temperature Sb surface. A significant focus will be placed on
presenting a symmetry argument that prohibits the existence of the
conventional uniform charge-4e phase. Furthermore, a comparison of various
vestigial phases will be conducted, emphasizing their distinct STM signatures.
Additionally, the role of CDW disorder in stabilizing vestigial
superconducting phases will be discussed, particularly in the context of
commensurate systems. The necessity of CDW disorder for the stabilization of
these vestigial phases will be explored and elucidated.

##  II Symmetry Analysis

In this study, our focus is on the low-temperature PDW state observed on the
Sb-cleaved surface of the Kagome superconductor CVS. On this surface, multiple
CDW states have been reported above the superconducting phase. The first CDW
state, with a $2a_{0}$ periodicity, is observed at 94K Ortiz et al. (2020);
Zhao et al. (2021a); Jiang et al. (2021); Li et al. (2021). Its wavevectors,
shown as blue dots in Fig.1, preserve the 6-fold rotational symmetry. The
second CDW state, with a $4a_{0}$ periodicity in the XY-plane, emerges around
60K Zhao et al. (2021a); Chen et al. (2021a). Its wavevectors, shown as red
dots, break the 6-fold rotational symmetry. Both CDW ordering temperatures
significantly exceed the superconducting transition temperature $T_{c}\approx
2.3K$.

Figure 1: The Bragg peaks from the Kagome lattice structure (black),
$2a_{0}$-CDW (blue) and $4a_{0}$-CDW (red). Inside the superconducting phase,
four additional peaks (black circles at $\bf\pm Q_{1,2}$) are observed in STM.

Deep in the SC phase at 300mK, four additional CDW $\rho_{\bf Q_{i}}$ with
peaks at $\bf\pm Q_{1,2}$ have been observed in STM Chen et al. (2021a), which
are the signature for PDW $\Delta_{\bf Q_{i}}$. Given the existence of a
uniform SC component $\Delta_{0}$, these CDW wave vectors $\bf Q_{i}$ are the
same as the PDW wave vectors through $\rho_{\bf
Q_{i}}=\Delta_{0}^{*}\Delta_{\bf Q_{i}}$. The same wave vectors are observed
in the spatial variation of the SC gap magnitudeChen et al. (2021a). The exact
critical temperature for the low-temperature CDW is currently unknown.
However, experimental evidence suggests that it is below 4.2K Chen et al.
(2021a). The peak at $\bf\pm Q_{3}$ is a higher harmonics of the
unidirectional CDW, which exists even above the SC phaseChen et al. (2021a),
so it does not break additional translational symmetry.

Figure 2: The unfolded Brillouin (black line), folded Brillouin zone due to
$2a_{0}$-CDW (blue line) and folded Brillouin zone due to both $2a_{0}$ and
$4a_{0}$-CDW (red line). In the last folded Brillouin zone, the PDW (black
circle) is located at the corner.

When considering the low-temperature phases ($T<4K$) in the presence of the
$2a_{0}$ and $4a_{0}$ CDW, it is appropriate to treat them as static
backgrounds. This leads to a larger unit cell, requiring the Brillouin zone
(BZ) folding. In Fig.2, the original Brillouin zone is depicted as a solid
black hexagon. Upon folding the Brillouin zone according to the $2a_{0}$ CDW
(blue dots), the resulting Brillouin zone (indicated by the blue line) remains
a hexagon, but with a halved size. The PDWs (black circles) are now folded
onto the M points of the resulting Brillouin zone. Consequently, the PDWs
$\Delta_{\bf Q_{i}}$ and $\Delta_{\bf-Q_{i}}$ become physically identical,
breaking the same symmetry.

We now proceed to fold the BZ according to the unidirectional $4a_{0}$-CDW
(red dots). The resulting folded Brillouin zone is depicted as a red
rectangle. Interestingly, all four PDWs (originally $\bf\pm Q_{1,2}$) are now
located at the corners of the rectangle (i.e. the $(\pi,\pi)$ point). This
implies that these PDWs are physically equivalent, exhibiting the same
symmetry breaking.

To summarize, the low-temperature phase observed on the Sb-layer of CVS
features a non-degenerate PDW $\Delta_{{\bf Q}=(\pi,\pi)}$ and a uniform
superconducting order parameter $\Delta_{\bf 0}$. This differs from the well-
known ‘LO’-like phase Larkin and Ovchinnikov (1972), which exhibits a doubly-
degenerate PDW $\Delta_{\bf\pm Q}$. Furthermore, the low-temperature phase in
CVS is distinct from the ‘FF’-like phase Fulde and Ferrell (1964) as the
wavevector $\bf Q$ preserves time-reversal symmetry.

Figure 3: (Left) The conventional phase diagram hosting a vestigial charge-4e
phase from doubly degenerate PDWs $\Delta_{\pm\bf Q}$Berg et al. (2009). The
blue line denotes the superconducting transition while the red line denotes
the CDW transition. The remaining symmetry in each phase is included.
Depending on the commensurability of the CDW, CDW transition breaks U(1)
(incommensurate) or $\mathbb{Z}_{m}$ (commensurate) symmetry. (Right) Unlikely
phase diagram hosting a charge-4e phase using only a non-degenerate PDW. It
cannot be achieved since the discrete symmetry breaking (black line)
generically happens above the continuous symmetry breaking (blue line).

The non-degenerate PDW $\Delta_{\bf Q}$ introduces an interesting scenario for
vestigial phases. Traditionally, vestigial phases require two independent
superconducting order parameters. In the case of the doubly degenerate PDW
$\Delta_{\bf\pm Q}$, the vestigial charge-4e phase arises from the uniform
order parameter $\Delta^{4e}_{0}=\Delta_{\bf Q}\Delta_{\bf-Q}$ Berg et al.
(2009). In this phase, $\Delta^{4e}_{0}$ is ordered while the corresponding
CDW $\rho_{2{\bf Q}}=\Delta_{\bf Q}\Delta_{\bf-Q}^{*}$ is disordered. The
phase diagram is depicted in the left panel of Fig.3, where the blue line
represents the superconducting transition and the red line represents the CDW
transition. The remaining symmetries in each phase are also indicated. The
superconducting transition breaks a U(1) symmetry, while the CDW transition
can either break an additional U(1) symmetry if the CDW $\rho_{2{\bf Q}}$ is
incommensurate or a $\mathbb{Z}_{m}$ symmetry if the periodicity of the CDW is
m times the lattice spacing. To support the charge-4e phase, it is necessary
for the superconducting transition to occur at a higher critical temperature
than the CDW transition. In the case of an incommensurate CDW, this
requirement can be met in a clean system as long as the superconducting phase
has a larger stiffness Berg et al. (2009). However, for a commensurate CDW,
discrete symmetry breaking generally occurs at a higher $T_{c}$ than
continuous symmetry breaking in a clean system. In such cases, CDW disorder
can be required to suppress the CDW $T_{c}$, allowing the charge-4e phase to
persist.

In the case of CVS, the conventional construction of the charge-4e phase is
not applicable due to the presence of the $2a_{0}$-CDW, which makes
$\Delta_{\bf\pm Q}$ the same order parameter. In principle, one could consider
a charge-4e phase where $\Delta_{\bf Q}^{2}$ is ordered while $\Delta_{\bf Q}$
remains disordered. The corresponding phase diagram is depicted in the right
panel of Fig.3. However, this phase cannot be stabilized. In this charge-4e
phase, $\Delta_{\bf Q}$ becomes an Ising variable, leading to the breaking of
the continuous U(1) symmetry into a $\mathbb{Z}_{2}$ symmetry. For this
charge-4e phase to be realized, it would require the U(1) continuous symmetry
breaking to occur at a higher temperature than the $\mathbb{Z}_{2}$ discrete
symmetry breaking for the same order parameter. However, in practice, this is
unlikely to happen since discrete symmetry breaking, which is less affected by
fluctuations, typically occurs at higher temperatures than continuous symmetry
breaking.

Vestigial phases, if existed, should have alternative order parameters,
constructed from two distinct SC orders. In this system, the choice is the PDW
$\Delta_{\bf Q}$ and the uniform SC $\Delta_{0}$. The vestigial charge-4e
order and the vestigial CDW order thus should be $\Delta^{4e}_{\bf
Q}=\Delta_{0}\Delta_{\bf Q}$ and $\rho_{\bf Q}=\Delta_{0}\Delta^{*}_{\bf Q}$.

As shown in the left panel of Fig.4, the phase diagram hosting the charge-4e
phase has the same essence as the conventional phase diagram in Fig.3. The
charge-4e SC transition (blue line) is a U(1) symmetry breaking, where
$\Delta^{4e}_{\bf Q}=\Delta_{0}\Delta_{\bf Q}$ is ordered. The CDW transition
(red line) is a $\mathbb{Z}_{2}$ (since ${2\bf Q}={\bf 0}$) translational
symmetry breaking, where $\rho_{\bf Q}=\Delta_{\bf Q}\Delta_{0}^{*}$ is
ordered. It is worth emphasizing that, the CDW here and below, corresponds to
peaks at $\pm Q_{1,2}$. Its critical temperature is below 4.2K. To host the
charge-4e phase, the continuous U(1) symmetry breaking needs to have a higher
$T_{c}$ than the discrete $\mathbb{Z}_{2}$ symmetry breaking. This is
difficult to achieve for a clean system. However, CDW disorder can suppress
the CDW $T_{c}$, allowing the charge-4e phase to survive. A numerical analysis
of the effect of the disorder can be found in Sec.IV.

Figure 4: Possible phase diagrams with the uniform SC $\Delta_{0}$ and the
non-degenerate PDW $\Delta_{\bf Q}$ at low temperature. The blue line is a
U(1) superconducting transition while the red line is a $\mathbb{Z}_{2}$ CDW
transition. (Left) The vestigial charge-4e SC phase is developed first, with
$\Delta_{0}$ and $\Delta_{\bf Q}$ coexisting in the charge-2e SC phase.
(Middle) The vestigial CDW phase is developed first, with $\Delta_{0}$ and
$\Delta_{\bf Q}$ coexisting in the charge-2e SC phase. (Right) No vestigial
phases are developed, with $\Delta_{0}$ and $\Delta_{\bf Q}$ ordered at
different temperatures.

##  III Experimental signatures

From the symmetry argument in the previous section, vestigial phases need to
have a composite order parameter from both $\Delta_{0}$ and $\Delta_{\bf Q}$.
In this section, we will illustrate the experimental signatures of possible
phase diagrams.

We will focus on the STM probe. In the charge-4e phase, there is a non-zero
density of states for gapless excitations, due to the lack of $\bf Q=0$
pairing components. This has the essence as the previous studies on the
doubly-degenerate PDW phaseLee (2014). As usual, the STM gap opening, i.e.
vanishing density of states at the Fermi level, is a signature for the
ordering of $\Delta_{0}$. And the CDW peaks at $\pm\bf Q_{1,2}$ are signatures
for the ordering of $\rho_{\bf Q}=\Delta_{0}\Delta_{\bf Q}^{*}$. We will use
these two signatures to distinguish phase diagrams with/without vestigial
phases.

In the vestigial SC phase $\Delta^{4e}_{\bf Q}\neq 0$ (left panel of Fig.4),
the PDW and the uniform charge-2e order can mutually induce each other through
the term $\Delta^{4e}_{\bf Q}\Delta^{*}_{0}\Delta^{*}_{\bf Q}$ in Landau
theory. Consequently, the gap opening for $\Delta_{\bf 0}$ and $\Delta_{\bf
Q}$ must happen at the same temperature. The STM gap opening and the CDW peaks
at $\pm\bf Q_{1,2}$ will thus appear at the same temperature (red line).

If the vestigial charge-4e phase does not exist, there could be a vestigial
CDW phase (middle panel). In this phase, $\rho_{\bf Q}=\Delta_{\bf
Q}\Delta_{\bf 0}^{*}$ is ordered, while SC is not. At the SC $T_{c}$, both
gaps $\Delta_{\bf 0}$ and $\Delta_{\bf Q}$ open for the same reason above. In
STM, the CDW peaks $\pm\bf Q_{1,2}$ (red line) will thus appear at a higher
temperature than the SC gap opening (blue line).

If neither vestigial phases exist, then the gaps $\Delta_{\bf 0}$ and
$\Delta_{\bf Q}$ will generically open at different temperatures, as shown in
the right panel of Fig.4. We take $\Delta_{0}$ as a stronger ordering with a
higher critical temperature as $|\Delta_{0}|$ is comparable to the
superconducting $T_{c}$ Chen et al. (2021a). In STM, SC gap opening (blue
line) will appear at a higher temperature than CDW peaks $\pm\bf Q_{1,2}$ (red
line).

For simplicity, we only discussed the above three possibilities. Other
possibilities (e.g. charge-6e phase) will lead to more complicated phase
diagrams. As long as vestigial SC phases (4e or 6e) exist, $\Delta_{0}$ and
$\Delta_{\bf Q}$ can induce each other in that phase. So the above properties
(same $T_{c}$ for gap opening and CDW peaks) still provide a straightforward
way to check the possibility of vestigial SC phases.

Another more direct way to detect the charge-4e phase is by phase-sensitive
measurement. In real space, the charge-4e order parameter reads:
$\Delta^{4e}({\bf r})=|\Delta^{4e}_{\bf Q}|\exp(i\bf Q\cdot r)$. When putting
this charge-4e SC next to a uniform charge-2e superconductor $\Delta^{2e}$, it
will induce a PDW order $\Delta_{\bf Q}=\Delta^{4e}_{\bf Q}\Delta^{2e,*}$, as
well as a CDW order $\rho_{\bf Q}=\Delta^{4e}_{\bf
Q}\Delta^{2e,*}\Delta^{2e,*}$ at the boundary.

##  IV  Disorder Effect

For the vestigial SC state to appear, the critical temperature $T_{SC}$ needs
to be higher than $T_{CDW}$. A commensurate CDW breaks discrete symmetry while
SC breaks continuous symmetry. In a clean system, thermal fluctuation in SC
would likely be stronger than CDW, leading to $T_{SC}<T_{CDW}$. In this case,
CDW disorder can help suppress $T_{CDW}$, thus allowing the vestigial SC phase
to survive. The effect can be described by the following Landau-Ginzburg-
Wilson free energy density:

$\begin{split}f({\bf
r})&=k_{\Delta}|\nabla\Delta|^{2}+\frac{u_{\Delta}}{2}|\Delta|^{2}+\frac{1}{4}|\Delta|^{4}\\\
&+k_{\phi}(\nabla\phi)^{2}+\frac{u_{\phi}}{2}\phi^{2}+\frac{1}{4}\phi^{4}+\frac{\lambda}{2}|\Delta|^{2}\phi^{2}+h({\bf
r})\phi\end{split}$ (1)

Here $\Delta$ is a complex field, describing the U(1) symmetry breaking for
the charge-4e SC. $\phi$ is a real field, describing the $\mathbb{Z}_{2}$
symmetry breaking for the commensurate CDW near the superconducting
transition. The $\lambda$ term describes the competition between them.

We include the CDW disorder as a real random field variable $h({\bf r})$, with
the spatially independent $h({\bf r})$ uniformly distributed within $[-D,D]$.
The parameter D then measures the disorder strength. Note that the SC cannot
be coupled to any random-field disorder due to the gauge symmetry.

This minimal model captures the most straightforward consequence of CDW
disorder, which is suppressing $T_{CDW}$ and thus allowing the vestigial SC
phase. We find it unnecessary to include all symmetry-allowed terms,
originated from the rotational symmetry breaking and the possible time-
reversal symmetry breaking in the normal phase. For example, the inclusion of
time-reversal symmetry breaking, represented as an internal magnetic field,
can lead to the formation of magnetic vortices within the superconducting
state. The introduction of random-$T_{c}$ disorder on the superconductivity
can also result in the presence of local SC patches, which can enhance the
superconducting $T_{c}$ under proper conditions. It is essential to clarify
that our intention in this work is not to undermine the importance of these
physics. Rather, our primary objective is to emphasize that, CDW disorder can
help the existence of the vestigial SC phase through its most straightforward
and non-exotic mechanism.

The random-field disorder is known to destroy long-range ordering in 2D, but
more realistically, we will focus on a finite-sized system. We compute the
correlation functions:
$c_{\Delta}(r)=\frac{1}{L^{2}}\sum_{x,y}\langle\Delta(x,y)\Delta^{*}(x+r,y+r)\rangle$
and
$c_{\phi}(r)=\frac{1}{L^{2}}\sum_{x,y}\langle\phi(x,y)\phi(x+r,y+r)\rangle$.
Here $\langle\rangle$ denotes the thermal average. Then we obtain
configurational average $\overline{c_{\Delta,\phi}(r)}$ over different
disorder samplings.

For a $L\times L$ periodic system, we mainly look at the correlation functions
at the longest distance
$C_{\Delta,\phi}\equiv\overline{c_{\Delta,\phi}(L/2)}$. The phase boundaries
can be determined when $C_{\Delta,\phi}$ reaches a threshold value.
$u_{\Delta}=-1$, $u_{\phi}=-0.95$ are taken for the classical Monte Carlo
calculation. Details can be found in the appendix.

Figure 5: Correlation functions at the longest distance, for (left) charge-4e
SC and (right) CDW, as a function of disorder strength $D$ and temperature
$T$.

Let us start with zero disorder $D=0$. At zero temperature (without thermal
fluctuation), SC should be a stronger order as we choose
$|u_{\Delta}|>|u_{\phi}|$. With thermal fluctuation (Fig.5), the commensurate
CDW has a higher $T_{c}$, since the discrete symmetry breaking suffers from
less fluctuation. This weak disorder limit corresponds to the phase diagram in
the middle panel of Fig.4, which has a vestigial CDW phase.

As a strong disorder suppresses the CDW, the charge-4e SC has a higher
$T_{c}$. This strong disorder limit corresponds to the left panel of Fig.4,
with a vestigial charge-4e phase. There is a tetracritical point (crossing
point between the two continuous phase transitions) in the disorder-
temperature phase diagram, separating the above two limits. In sum, CDW
disorder can be crucial for stabilizing the vestigial charge-4e phase in a
general commensurate system.

##  V Discussion

In this study, our focus is specifically on the low-temperature surface
properties of CsV3Sb5. Our analysis is solely based on the existing scanning
tunneling microscopy results without making any additional assumptions on the
3D wave vectors of the CDW and PDW. On the surface of CVS, the PDW state is
found to be non-degenerate, implying that there is only one unique approach to
constructing a vestigial charge-4e phase. By utilizing the experimental STM
data, we aim to provide a comprehensive understanding of the surface behavior
and the specific characteristics of the vestigial charge-4e phase in CVS.

A natural continuation of this study is to investigate the behavior of CsV3Sb5
in bulk. It is crucial to explore the existence of a three-dimensional PDW
state and determine its out-of-plane wavevector to understand the symmetry-
breaking phenomena. For instance, one may question whether the points
$\pm\mathbf{Q}_{i}$ remain equivalent in the three-dimensional case.

Observations have shown that the $4a_{0}$-CDW is confined to the Sb surface
Zhao et al. (2021a); Chen et al. (2021a); Liang et al. (2021), indicating that
it is unlikely to be a bulk effect. Consequently, the PDW states at
wavevectors $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ are expected to differ. In
the bulk, with the presence of more than two superconducting order parameters,
multiple scenarios can give rise to vestigial charge-4e phases. By
constructing the corresponding surface phase diagram based on a given bulk
phase diagram, STM probes can play a crucial role in ruling out certain
possibilities in the bulk. A specific example illustrating this process can be
found in the appendix of this paper. However, it is important to note that our
analysis of the bulk PDW wavevectors involves certain assumptions, as there is
currently a lack of experimental results regarding the 3D wave vectors of the
PDW in CsV3Sb5.

Although this study primarily focuses on the Kagome superconductor CsV3Sb5,
similar investigations could be extended to the other two Kagome
superconductors, KV3Sb5 and RbV3Sb5. In KV3Sb5 and RbV3Sb5, the $2a_{0}$ CDW
has been observed through various experimental techniques Yin et al. (2021);
Ortiz et al. (2021b); Jiang et al. (2021); Li et al. (2022b, 2021); Miao et
al. (2021); Liu et al. (2021); Luo et al. (2022b). The unidirectional $4a_{0}$
CDW observed in CsV3Sb5 is not detected in these two materials. To date, no
PDW has been observed in KV3Sb5 and RbV3Sb5, but there have been reports of
dispersive scattering wave vectors at $\bf Q_{i}$ in the normal state of
KV3Sb5 Li et al. (2023).

To fully understand the superconducting states in these two materials, it is
thus crucial to check the existence of any non-dispersive CDW peaks at low
temperature. For instance, if the dispersive scattering in KV3Sb5 became non-
dispersive at $T_{CDW}$, this would suggest the existence of PDW in KVS.
Interestingly, due to the presence of the $2a_{0}$ CDW in that system, CDW
peaks at $\pm Q_{1}$ are equivalent. So the conventional uniform vestigial
superconductivity $\Delta_{Q_{1}}\Delta_{-Q_{1}}$ should not exist. The
correct form of the vestigial superconductivity thus depends on results from
future experiments.

Acknowledgments - We thank S. Raghu, S. Kivelson, D. Agterberg, and Z. Han for
helpful discussions.

## References

  * Ortiz et al. (2019) B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, et al., Physical Review Materials 3, 094407 (2019).
  * Ortiz et al. (2020) B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, et al., Phys. Rev. Lett. 125, 247002 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.125.247002.
  * Jiang et al. (2023) K. Jiang, T. Wu, J.-X. Yin, Z. Wang, M. Z. Hasan, S. D. Wilson, X. Chen, and J. Hu, National Science Review 10, nwac199 (2023).
  * Xu et al. (2021) H.-S. Xu, Y.-J. Yan, R. Yin, W. Xia, S. Fang, Z. Chen, Y. Li, W. Yang, Y. Guo, and D.-L. Feng, Physical Review Letters 127, 187004 (2021).
  * Chen et al. (2021a) H. Chen, H. Yang, B. Hu, Z. Zhao, J. Yuan, Y. Xing, G. Qian, Z. Huang, G. Li, Y. Ye, et al., Nature 599, 222 (2021a).
  * Liang et al. (2021) Z. Liang, X. Hou, F. Zhang, W. Ma, P. Wu, Z. Zhang, F. Yu, J.-J. Ying, K. Jiang, L. Shan, et al., Physical Review X 11, 031026 (2021).
  * Wang et al. (2021a) Z. Wang, Y.-X. Jiang, J.-X. Yin, Y. Li, G.-Y. Wang, H.-L. Huang, S. Shao, J. Liu, P. Zhu, N. Shumiya, et al., Physical Review B 104, 075148 (2021a).
  * Li et al. (2021) H. Li, T. T. Zhang, T. Yilmaz, Y. Y. Pai, C. E. Marvinney, A. Said, Q. W. Yin, C. S. Gong, Z. J. Tu, E. Vescovo, et al., Phys. Rev. X 11, 031050 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.031050.
  * Miao et al. (2021) H. Miao, H. X. Li, W. R. Meier, A. Huon, H. N. Lee, A. Said, H. C. Lei, B. R. Ortiz, S. D. Wilson, J. X. Yin, et al., Phys. Rev. B 104, 195132 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.104.195132.
  * Ortiz et al. (2021a) B. R. Ortiz, S. M. L. Teicher, L. Kautzsch, P. M. Sarte, N. Ratcliff, J. Harter, J. P. C. Ruff, R. Seshadri, and S. D. Wilson, Phys. Rev. X 11, 041030 (2021a), URL https://link.aps.org/doi/10.1103/PhysRevX.11.041030.
  * Song et al. (2022) D. Song, L. Zheng, F. Yu, J. Li, L. Nie, M. Shan, D. Zhao, S. Li, B. Kang, Z. Wu, et al., Science China Physics, Mechanics & Astronomy 65, 247462 (2022).
  * Jiang et al. (2021) Y.-X. Jiang, J.-X. Yin, M. M. Denner, N. Shumiya, B. R. Ortiz, G. Xu, Z. Guguchia, J. He, M. S. Hossain, X. Liu, et al., Nature Materials 20, 1353 (2021).
  * Zhao et al. (2021a) H. Zhao, H. Li, B. R. Ortiz, S. M. Teicher, T. Park, M. Ye, Z. Wang, L. Balents, S. D. Wilson, and I. Zeljkovic, Nature 599, 216 (2021a).
  * Mielke III et al. (2022) C. Mielke III, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. Lei, J. Chang, et al., Nature 602, 245 (2022).
  * Yu et al. (2021a) L. Yu, C. Wang, Y. Zhang, M. Sander, S. Ni, Z. Lu, S. Ma, Z. Wang, Z. Zhao, H. Chen, et al., arXiv preprint arXiv:2107.10714 (2021a).
  * Li et al. (2022a) H. Li, S. Wan, H. Li, Q. Li, Q. Gu, H. Yang, Y. Li, Z. Wang, Y. Yao, and H.-H. Wen, Physical Review B 105, 045102 (2022a).
  * Feng et al. (2021a) X. Feng, K. Jiang, Z. Wang, and J. Hu, Science bulletin 66, 1384 (2021a).
  * Park et al. (2021) T. Park, M. Ye, and L. Balents, Physical Review B 104, 035142 (2021).
  * Feng et al. (2021b) X. Feng, Y. Zhang, K. Jiang, and J. Hu, Physical Review B 104, 165136 (2021b).
  * Denner et al. (2021) M. M. Denner, R. Thomale, and T. Neupert, Physical Review Letters 127, 217601 (2021).
  * Lin and Nandkishore (2021) Y.-P. Lin and R. M. Nandkishore, Physical Review B 104, 045122 (2021).
  * Yu et al. (2021b) F. Yu, T. Wu, Z. Wang, B. Lei, W. Zhuo, J. Ying, and X. Chen, Physical Review B 104, L041103 (2021b).
  * Xiang et al. (2021) Y. Xiang, Q. Li, Y. Li, W. Xie, H. Yang, Z. Wang, Y. Yao, and H.-H. Wen, Nature communications 12, 6727 (2021).
  * Ni et al. (2021) S. Ni, S. Ma, Y. Zhang, J. Yuan, H. Yang, Z. Lu, N. Wang, J. Sun, Z. Zhao, D. Li, et al., Chinese Physics Letters 38, 057403 (2021).
  * Ratcliff et al. (2021) N. Ratcliff, L. Hallett, B. R. Ortiz, S. D. Wilson, and J. W. Harter, Physical Review Materials 5, L111801 (2021).
  * Wang et al. (2021b) Z. Wang, Q. Wu, Q. Yin, C. Gong, Z. Tu, T. Lin, Q. Liu, L. Shi, S. Zhang, D. Wu, et al., Physical Review B 104, 165110 (2021b).
  * Luo et al. (2022a) J. Luo, Z. Zhao, Y. Zhou, J. Yang, A. Fang, H. Yang, H. Gao, R. Zhou, and G.-q. Zheng, npj Quantum Materials 7, 30 (2022a).
  * Li et al. (2023) H. Li, D. Oh, M. Kang, H. Zhao, B. R. Ortiz, Y. Oey, S. Fang, Z. Ren, C. Jozwiak, A. Bostwick, et al., arXiv preprint arXiv:2303.07254 (2023).
  * Mu et al. (2021) C. Mu, Q. Yin, Z. Tu, C. Gong, H. Lei, Z. Li, and J. Luo, Chinese Physics Letters 38, 077402 (2021).
  * Yu and Li (2012) S.-L. Yu and J.-X. Li, Physical Review B 85, 144402 (2012).
  * Kiesel et al. (2013) M. L. Kiesel, C. Platt, and R. Thomale, Physical review letters 110, 126405 (2013).
  * Wang et al. (2013) W.-S. Wang, Z.-Z. Li, Y.-Y. Xiang, and Q.-H. Wang, Physical Review B 87, 115135 (2013).
  * Duan et al. (2021) W. Duan, Z. Nie, S. Luo, F. Yu, B. R. Ortiz, L. Yin, H. Su, F. Du, A. Wang, Y. Chen, et al., Science China Physics, Mechanics & Astronomy 64, 107462 (2021).
  * Luo et al. (2022b) H. Luo, Q. Gao, H. Liu, Y. Gu, D. Wu, C. Yi, J. Jia, S. Wu, X. Luo, Y. Xu, et al., Nature communications 13, 273 (2022b).
  * Zhao et al. (2021b) C. Zhao, L. Wang, W. Xia, Q. Yin, J. Ni, Y. Huang, C. Tu, Z. Tao, Z. Tu, C. Gong, et al., arXiv preprint arXiv:2102.08356 (2021b).
  * Wu et al. (2021) X. Wu, T. Schwemmer, T. Müller, A. Consiglio, G. Sangiovanni, D. Di Sante, Y. Iqbal, W. Hanke, A. P. Schnyder, M. M. Denner, et al., Physical review letters 127, 177001 (2021).
  * Gu et al. (2022) Y. Gu, Y. Zhang, X. Feng, K. Jiang, and J. Hu, Physical Review B 105, L100502 (2022).
  * Gupta et al. (2022) R. Gupta, D. Das, C. H. Mielke III, Z. Guguchia, T. Shiroka, C. Baines, M. Bartkowiak, H. Luetkens, R. Khasanov, Q. Yin, et al., npj Quantum Materials 7, 49 (2022).
  * Chen et al. (2021b) K. Chen, N. Wang, Q. Yin, Y. Gu, K. Jiang, Z. Tu, C. Gong, Y. Uwatoko, J. Sun, H. Lei, et al., Physical Review Letters 126, 247001 (2021b).
  * Yu et al. (2021c) F. Yu, D. Ma, W. Zhuo, S. Liu, X. Wen, B. Lei, J. Ying, and X. Chen, Nature communications 12, 3645 (2021c).
  * Chen et al. (2021c) X. Chen, X. Zhan, X. Wang, J. Deng, X.-B. Liu, X. Chen, J.-G. Guo, and X. Chen, Chinese Physics Letters 38, 057402 (2021c).
  * Zhang et al. (2021) Z. Zhang, Z. Chen, Y. Zhou, Y. Yuan, S. Wang, J. Wang, H. Yang, C. An, L. Zhang, X. Zhu, et al., Physical Review B 103, 224513 (2021).
  * Berg et al. (2009) E. Berg, E. Fradkin, and S. A. Kivelson, Nature Physics 5, 830 (2009).
  * Ge et al. (2022) J. Ge, P. Wang, Y. Xing, Q. Yin, H. Lei, Z. Wang, and J. Wang, arXiv preprint arXiv:2201.10352 (2022).
  * Larkin and Ovchinnikov (1972) A. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP 34, 1144 (1972).
  * Fulde and Ferrell (1964) P. Fulde and R. A. Ferrell, Physical Review 135, A550 (1964).
  * Lee (2014) P. A. Lee, Physical Review X 4, 031017 (2014).
  * Yin et al. (2021) Q. Yin, Z. Tu, C. Gong, Y. Fu, S. Yan, and H. Lei, Chinese Physics Letters 38, 037403 (2021).
  * Ortiz et al. (2021b) B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. Teicher, R. Seshadri, and S. D. Wilson, Physical Review Materials 5, 034801 (2021b).
  * Li et al. (2022b) H. Li, H. Zhao, B. R. Ortiz, T. Park, M. Ye, L. Balents, Z. Wang, S. D. Wilson, and I. Zeljkovic, Nature Physics 18, 265 (2022b).
  * Liu et al. (2021) Z. Liu, N. Zhao, Q. Yin, C. Gong, Z. Tu, M. Li, W. Song, Z. Liu, D. Shen, Y. Huang, et al., Physical Review X 11, 041010 (2021).

##  Appendix A Methods

In this section, we show the detailed steps for the classical Monte Carlo
calculation.

We start with a $L\times L$ square lattice with periodic boundary conditions.
At each site $\bf r$, we initialize the CDW disorder $h({\bf r})$. The
disorder is chosen independently from a uniform distribution $[-D,D]$. We also
initialize a random configuration for $\Delta({\bf r})$ and $\phi({\bf r})$.
Here, $\Delta({\bf r})$ are complex variables while $\phi({\bf r})$ are real
variables.

We can now compute the free energy, for a given $\Delta({\bf r})$ and
$\phi({\bf r})$ configuration. The kinetic terms live on the bond. For
example, the kinetic term for $\phi$ is $\sum_{\langle{\bf
rr^{\prime}}\rangle}\frac{k_{\phi}}{2}(\phi({\bf r})-\phi({\bf
r^{\prime}}))^{2}$, summed up nearest neighbored sites $\langle{\bf
rr^{\prime}}\rangle$. The other terms live on sites. For example, the disorder
term is $\sum_{\bf r}h({\bf r})\phi({\bf r})$, summed up all sites.

We now propose a random change for $\phi({\bf r})$ at a particular site ${\bf
r}$. If the change in the free energy $\Delta E$ is negative, then this change
is accepted. Otherwise, its acceptance rate is $\exp(-\Delta E/T)$. And then
we move to another site. As we go through all $\phi({\bf r})$, we will perform
the same consideration for the complex variables $\Delta({\bf r})$. As we go
through all $\phi({\bf r})$ and $\Delta({\bf r})$, we finish one ‘step’ in the
simulation.

After a sufficient number of steps, the system is stabilized and we are ready
to perform measurements, as written in the main text. Two measurements need to
be separated by a sufficient number of steps, to make them independent. By
averaging out these measurements, we obtain the thermal average
$\langle...\rangle$.

In the calculations, we take $L=20$, $k_{\Delta}=k_{\phi}=1$, $u_{\Delta}=-1$,
$u_{\phi}=-0.95$ and $\lambda=0.2$. Classical Monte Carlo calculation is
performed with $10^{6}$ steps and $10^{3}$ measurements.

We generate different 40 disorder configurations. By averaging their
correlation function, we obtain the final disorder average.

##  Appendix B An example for the bulk

In the main text, we have identified three distinct surface phase diagrams
that can be differentiated using STM probes. In this section, we outline the
methodology for deriving a surface phase diagram based on a given bulk phase
diagram.

Since the 3D wave vectors of the bulk PDW are unknown, we have to make some
assumptions about them. For simplicity, we assume six PDW
$\pm\mathbf{Q}_{1,2,3}$ with zero out-of-plane wave vectors. Due to the lack
of the $4a_{0}$ CDW in the bulk, $\mathbf{Q}_{1,2}$ are no longer equivalent
points. And $\mathbf{Q}_{3}$ is unequivalent to $\bf 0$.

Figure 6: An example of obtaining the surface phase diagram from the bulk
phase diagram.

Due to the $2a_{0}$ CDW (no matter whether it is $2\times 2\times 2$ or
$2\times 2\times 4$), $\pm\mathbf{Q}_{i}$ are equivalent points because of the
above assumptions. Compared with the surface, the bulk does not have the
$4a_{0}$ CDW. So the total symmetry is larger. It becomes
U(1)$\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}$. Here U(1) describes the
superconductivity. The first $\mathbb{Z}_{2}$ symmetry (denoted as
$\mathbb{Z}_{2}^{(a)}$ below) is the same as in the main text. The second
$\mathbb{Z}_{2}$ symmetry (denoted as $\mathbb{Z}_{2}^{(b)}$ below)
corresponds to the $4a_{0}$ CDW on the surface.

Now it is clear how to get a surface phase diagram from a bulk phase diagram.
One only needs to further break an additional $\mathbb{Z}_{2}^{(b)}$ symmetry.

In the following example, the three symmetries are assumed to be broken at
different temperatures. As we go to the surface phase diagram, the
$\mathbb{Z}_{2}^{(b)}$ symmetry is no longer present due to the $4a_{0}$ CDW.
The resulting surface phase diagram is the middle one in Fig.4 in the main
text.

In the presence of multiple PDW states in the bulk, there exists a variety of
possible bulk phase diagrams. However, by employing the STM probe technique
proposed in the main text, it becomes feasible to identify the correct surface
phase diagram. This, in turn, provides valuable insights into the potential
bulk phase diagrams.


# Rapid design of fully soft deployable structures via kirigami cuts and
active learning

Leixin Ma1 Mrunmayi Mungekar1 Vwani Roychowdhury2 V.R<EMAIL_ADDRESS>M.K.
Jawed1 M.K.J<EMAIL_ADDRESS>1Dept. of Mechanical and Aerospace
Engineering, University of California, Los Angeles, CA 90095, USA
2Dept. of Electrical and Computer Engineering, University of California, Los
Angeles, CA 90095, USA

###### Abstract

Soft deployable structures – unlike conventional piecewise rigid deployables
based on hinges and springs – can assume intricate 3-D shapes, thereby
enabling transformative technologies in soft robotics, shape-morphing
architecture, and pop-up manufacturing. Their virtually infinite degrees of
freedom allow precise control over the final shape. The same enabling high
dimensionality, however, poses a challenge for solving the inverse design
problem involving this class of structures: to achieve desired 3D structures
it typically requires manufacturing technologies with extensive local
actuation and control during fabrication, and a trial and error search over a
large design space. We address both of these shortcomings by first developing
a simplified planar fabrication approach that combines two ingredients: strain
mismatch between two layers of a composite shell and kirigami cuts that
relieves localized stress. In principle, it is possible to generate targeted
3-D shapes by designing the appropriate kirigami cuts and selecting the right
amount of prestretch, thereby eliminating the need for local control. Second,
we formulate a data-driven physics-guided framework that reduces the
dimensionality of the inverse design problem using autoencoders and
efficiently searches through the “latent” parameter space in an active
learning approach. We demonstrate the effectiveness of the rapid design
procedure via a range of target shapes, such as peanuts, pringles, flowers,
and pyramids. Tabletop experiments are conducted to fabricate the target
shapes. Experimental results and numerical predictions from our framework are
found to be in good agreement.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

Morphing planar geometry to 3D shapes can find a wide variety of engineering
applications [1] from additive and subtractive manufacturing to soft actuators
[2, 3] and architecture [4, 5, 6]. Multiple mechanisms have been reported for
the 2D to 3D transformation, including residual stress-induced bending [7],
temperature-induced growth [8], inflatable membranes [9], composite materials
controlled by external stimuli (e.g., temperature and pH) [1], paper folding
[10], swelling [11], mechanical loads, and boundary conditions [12, 13, 14].
The fabrication of these shape-morphing structures often requires detailed
local control of the geometry, curvature, and stress [15, 1, 9, 16, 17, 7, 13,
14]. Even though optimizing and realizing arbitrary 3D deformed shapes, such
as a human face, is possible, this is often done at a cost of complicated
fabrication technique [1]. Moreover, the inverse design problem of optimizing
the physical parameters to achieve the targeted shape typically requires a
trial and error search over a high-dimensional space. To address such
shortcomings, this work introduces a new paradigm for planar manufacturing,
where a soft kirigami composite (a bilayer shell) deforms from 2D plane to the
target 3D shape due to the coupling between two key mechanisms, i.e., kirigami
(i.e. material removal) [18, 3, 19] and strain mismatch [20, 21]. Compared to
fabrication requiring precise control of various parts of the structure, the
reported technique first developed in [22] – described later in this paper –
is much easier to implement.

The convenience of the proposed manufacturing technique, however, is not
enough to guarantee adoption unless the associated high-dimensional inverse
design problem can also be solved efficiently. In particular, we next address
how one can efficiently optimize the structural and geometrical properties to
achieve a target 3D shape. In order to define the problem let us first
consider the case without any kirigami cuts. It has been shown that by varying
only the prestretch (but no kirigami cuts), a planar composite plate with
strain mismatch can morph into different 3D shapes; these shapes, however, are
limited to only the few mode shapes [20, 22]. To achieve programmability of
more classes of shapes one can locally remove materials, thereby introducing
kirigami cuts. The optimization problem then reduces to finding kirigami
patterns in 2D (i.e. marking the areas where material has to be removed) that
would morph the planar structure to a desired 3D shape. However, kirigami cuts
have highly nonlocal impact on the structure. A cut at a certain location may
affect the global shape of the deformed shell. Hence, the kirigami pattern
optimization for target 3D shapes cannot be conducted locally in space, but
requires a global approach to explore the large design space, which is a
common problem in many inverse metamaterial designs [12, 23].

Such inverse design problems where metamaterials are designed for target
functionalities have received significant boost from recent advances in
machine learning (ML). Various ML algorithms, such as Variational Autoencoder
(VAE) [24, 25] and Generative Adversarial Network [26, 23] have been
successfully applied. However, most of these ML-aided inverse design methods
still depend on the generation of a computationally prohibitive number of
forward simulation data [27, 26]. A network is then trained to learn the
inverse map in a supervised manner: the output of the forward simulation is
used as input, and the corresponding design parameters as the desired output.
This computational overhead of generating forward simulation data is
particularly severe for our case. Given a kirigami pattern and other
parameters, it takes minutes of computing time to generate corresponding 3D
shapes. To train an inverse ML system, one could require millions of such
forward simulations, making such inverse approaches computationally
intractable for our application.

In this paper, we make a judicious use of recent advances in dimension
reduction techniques (e.g., VAEs), and integrate them with active learning
techniques (e.g., Bayesian optimization), where forward computations are done
on-demand; that is, only for local explorations. Such active learning
techniques, when used by themselves, typically do not work well for high
dimensional search spaces such as in our application. The use of model
reduction techniques of ML make these methods practicable. First, classes of
candidate kirigami patterns that possess the same symmetry property as the
target 3D shape are generated. A VAE is used to parameterize symmetric
patterns with low-dimensional latent features and, thus, reduce the dimension
of the inverse problem. Then, the optimal latent features are searched
iteratively using the Bayesian Optimization framework [28]. Numerical
simulations and controlled tabletop experiments are conducted to create target
3D shapes with different symmetric properties. Excellent agreement is found
between experiments and simulations, which demonstrates the prospect of using
this method in algorithmic design of metamaterials.

Figure 1: The inverse problem aims at designing target 3D shapes by finding
the optimal fabrication strategies. Six target shapes are shown on left, while
the six experimental results correspond to the target shapes shown on the
right. The planar fabrication uses a kirigami layer and a bottom layer is
shown in the middle. The three classes of kirigami layers correspond to 1)
unidirectional strips 2) reflectional symmetry 3) four-fold radial symmetry.
Figure 2: A. Overview of the fabrication concept: The composite structure with
a kirigami layer and a substrate layer bend to a free buckling shape under
radial prestretch. B. Flow chart of the data-driven design and optimization of
soft kirigami composite (B1) A VAE to reduce the dimension of kirigami
patterns to a small number of latent variables (B2) A Bayesian Optimization
loop that iteratively searches the optimal combination of latent kirigami
pattern variables, size of the structure, and amount of prestretch that
results to the target 3D topology.

## Overall concept

### Problem description

The inverse problem aims at designing target 3D shapes by finding the optimal
fabrication strategies, as shown in Fig. 1. The proposed fabrication technique
is shown in Fig. 2A1 to A4. We start from two thin flexible and stretchable
plates (Fig. 2A1). The maximum radius of the top layer is $R$. First, strain
mismatch between two plates or layers is created by radially stretching the
bottom one with the amount of prestretch $\lambda$ (Fig. 2A2). Then, the top
layer of the same radius as the stretched bottom layer is glued onto the
bottom layer (Fig. 2A3). The strain mismatch between the two layers induces
out-of-plane buckling (Fig. 2A4), since bending is less energetically
expensive than compression for thin shells [11]. However, using strain
mismatch in composite shells yields a limited number of structural modes [20,
22], and the elastic shells may easily wrinkle (localized deformation) to
relax the compressive stresses instead of a global change in shape [29].
Hence, material needs to be strategically removed from the top layer to create
kirigami patterns in order to yield a target 3D shape. The theoretical
characterization of these kirigami-aided stress relief is still not fully
available. The material removal expands the possible number of attainable 3D
shapes from 2D, with the out-of-plane buckling depending on both the magnitude
of prestretch and geometric parameters (including kirigami cuts) [30, 31].
Based on this rapid fabrication technique, we aim at inversely designing the
optimal kirigami patterns, size of the structure (measured by radius $R$), and
prestretch $\lambda$ such that a target 3D shape can be achieved.

### Machine learning framework for design of soft kirigami composite structure

To tackle the challenges of Edisonian searches, our machine learning-aided
design process is illustrated in Fig. 2B1 and B2, and discussed in more detail
below.

### Dimension reduction: Compact and continuous representation of kirigami
patterns

A kirigami pattern can be represented by an $N\times N$ binary image, where
the pixels in the uncut areas are represented by 1’s. This would correspond to
a very high ($N^{2}$) dimensional search space. Recent advances in Computer
Vision (CV), however, have shown that such a family of images typically lies
on a much smaller dimensional (say $D$) manifold, referred to as the latent
space. One can then search only over a $D$-dimensional space (e.g., $D$=6 in
our design), leading to considerable computational savings. A VAE is one such
computational model that can be used to learn a generative model of kiragami
patterns, where (a) every kirigami image can be mapped to a $D$-dimenional
vector using an Encoder network, and (b) every $D$-dimensional vector in the
latent space can be mapped to a sample kirigami image using a Decoder network.
The Encoder and Decoder networks are trained simultaneously, such that the
latent features are representative enough for various kirigami patterns. Once
trained, any active learning framework can view the space of kirigami patterns
as a $D$-dimensional continuous space. This allows one to train the VAE using
only a limited set of representative kirigami patterns, and using the VAE to
interpolate and generate potentially infinite number of similar kirigami
patterns. The detailed architecture and hyperparameters for the VAE are
described in the supplementary material.

### Design of numerical experiments

Optimal design parameters, including the kirigami latent features
($D$-dimensional), size of the structure (measured by radius), and prestretch
that yield the target 3D shape were searched by performing iterative Bayesian
Optimizations. To search global optimums, the optimization aims to “explore”
(prioritizing regions with large uncertainty) and “exploit” (focusing on the
regions with minimum loss function) the high dimensional parametric spaces,
until a suitable loss function is minimized [32].

To initialize the optimization process, we randomly sampled ten combinations
of latent features, sizes of the structure, and prestretches. For each
combination of the proposed kirigami patterns and prestretch, numerical
experiments were performed via finite element simulations. A Gaussian process
regression model was constructed to approximate the unknown effect of the
design variables on the loss function, which is the error of the 3D shape
between the simulations and design target. A common way to represent 3D data
is via the projected image of the height [33]. The negative structural
similarity index is chosen as the loss function, which is used to characterize
the dissimilarity between two images. Compared to mean squared error, the
structural similarity index (SSIM) is a perceptually-motivated loss function,
which is found to have better performance for image restoration tasks than the
squared $l_{2}$ norm of the error [34].

Once the Gaussian process model is constructed, an expected improvement
function [28], that increases with both the mean of the loss function and its
uncertainty, is calculated. Then, the combination of latent features, size of
the structure, and prestretch that maximizes the expected improvement function
is selected. The corresponding kirigami patterns for that combination of
latent features can be recovered from the trained decoder network. The
recovered patterns are used to create the mesh for the kirigami layer to
perform next round of finite element simulation. The negative SSIM between the
results of finite element simulation and target shape is used to update the
Gaussian process model and reassess the expected improvement function. Such a
process continues iteratively. Once the optimization converges, we verify its
effectiveness using precision desktop experiments.

Figure 3: A. Examples of different classes of kirigami patterns: 1)
unidirectional fibers, 2) shapes of reflectional symmetry with 2 lines of
symmetry, 3) four-fold radial symmetry. B. Experimental setup. (B1) Schematic
representation of the system consist of two-knob stages (1) on 250 mm tracks
(2) knobs to adjust the stage location (B2) Snapshot of the system when the
substrate layer (3) is stretched (B3) Attach a kirigami layer (4) on top of
the substrate layer (B4) Release the prestretch by cutting along the outline
of the circular substrate.

### Generation of candidate kirigami patterns for VAE training

We tested the performance of three classes of candidate kirigami patterns, as
shown in Fig. 3A1 to A3. The first class, as shown in Fig. 3A1, is made of
uni-directional fibers. For example, Hanakata et al. [24] demonstrated that
the uni-directional patterns have good expressive power, and can be used to
interpolate a variety of mixed kirigami cuts. The other classes of kirigami
patterns exhibit the same type of symmetric property as the target 3D shapes,
as shown in Fig. 3A2 and A3. We will see later in the paper (Fig. 5A and B)
that the peanut and pringle like shapes have reflectional symmetry with
respect to the two principal axes, while the flower shape has four-fold radial
symmetry in Fig. 5D. To create symmetric kirigami patterns, the kirigami cuts
are created in certain regions, with the patterns in the remaining regions
directly created via rotations or reflections (see _SI Text_). Both types of
non-symmetric and symmetric patterns are augmented via rotation. After
rotating the kirigami patterns two to four times, the total number of images
for each class is 40955, 40955, and 24573, respectively. These images are used
to train VAE models to find low-dimensional representations.

## Physical and numerical Experiments

### Desktop-scale physical experiments

Fig. 3B1 to B4 present our experimental setup and the key steps in
fabrication. The experimental setup consists of four linear translation stages
(250 mm travel, Thorlabs). The substrate and the kirigami layers are made of
hyper-elastic materials (see Materials & Methods). A large substrate layer is
applied with the four corners fixed to the four stages (Fig. 3B1). Then, to
impose the radial strain upon the substrate, we stretch the substrate layer
with the same amount in the horizontal and vertical directions (Fig. 3B2).
Once stretching is complete, a kirigami layer is glued on top of the substrate
(Fig. 3B3). We ensure that the central region of the substrate layer – onto
which the kirigami layer is glued – is uniformly stretched, not affected by
the boundaries. Afterwards, Fig. 3B4 shows that the excess substrate is cut
away along the outline of the circular substrate. The composite structure then
morphs to a certain 3D shape due to the strain mismatch in the two layers.

### Finite element-based numerical simulations

The finite element software, Abaqus [35], is used to model the nonlinear large
deformation of the hyper-elastic composite structures. In the simulations, we
use Mooney Rivlin model to simulate the hyper-elastic materials of the
substrate and the kirigami layer being tested. Each layer of the composite
structure was modeled as a shell to reduce the computational cost. The two
shell structures are constrained so that their normals match. The center of
the composite structure is fixed. The prestretch is modeled as isotropic
thermal expansion (see SI Text). The equivalent temperature corresponding to
the prestretch is computed. Then, the dynamics of the composite structure (for
the given candidate kirigami pattern) is modeled by Abaqus as the substrate is
cooled back to the reference temperature. The resulting 3D shape is obtained
after each such simulation converges. A nonlinear quasi-static simulation is
conducted for each design inputs.

The final 3D shape when the active learning search converges will be referred
to as the optimal predicted shape; the corresponding kirigami patterns,
radius, and prestretch are the predicted optimal parameters. The detailed
structural properties of the composite structure are included in the
“Materials and Methods” section. The numerical simulation is the most time-
consuming step, which takes about 1-3 minutes wall clock time on a desktop
computer (Ryzen 2950wx CPU @ 2.4 GHz). In total, 100 simulations take about 6
hours. Since forward simulations are expensive, cutting down the number of
runs is important. In this study, we are able to reduce the total number of
forward runs from millions to just around 100, while obtaining accurate target
shapes. Such computational tractability without compromising performance is
one of the key advantages of the proposed machine learning-aided framework.

## Results

In this section, we test drive the framework with a few example target shapes
with different distributions of curvatures and symmetric properties. During
the optimization, the range of radius of the kirigami structure is constrained
between 24 mm and 39 mm while the amount of prestretch $\lambda=1+\epsilon$ is
searched between 1.05 to 1.4, which means that the bottom layer is stretched
by 5% to 40% along both directions. Since in this range of prestretch, the
structural property is still properly modeled.

###  .1 Target 3D shapes with reflectional symmetry

As an exemplar of this class, we pick a bilobe structure, resembling a peanut
as the target shape. Automatically generated 40,955 (64 by 64) binary kirigami
images with reflectional symmetry are used to train the VAE. We find that just
a six dimensional latent space is able to reconstruct the training kirigami
data with very high accuracy (measured by SSIM $\approx$ 0.97). Thus, the
search space of kirigami patterns have been reduced from a $4098$-dimensional
binary image space to only a six dimensional continuous search space.New
kirigami patterns can be generated by the VAE, and they are also found to
exhibit reflectional symmetry, as demonstrated in the supplementary material.
More importantly, as demonstrated in Fig. 4 D, the interpolated new kirigami
patterns are often needed to generate 3D shapes that best match the target
shapes. We next illustrate some of the salient characteristics of our inverse
design framework.

Figure 4: A. Search trajectory Variation of SSIM between target peanut shape
and predicted free buckling shape over iterations. B. Multiple optimal
solutions B1-B4.The top four optimal design solutions, including the kirigami
patterns, prestretch and krigami size. The corresponding height distribution
for each design is presented at the bottom. C. Exploration during optimization
C1-C2.Two examples showing the algorithm explores new kirigami patterns, that
does not lead to closer solutions. D. Interpolated kirigami patterns are often
required to better approximate target 3D shapes:  D1-D2. Two nearest neighbors
for kirigami patterns found at iteration No. 78, whose free buckling shape
gives lower SSIM. E. Effect of initialization E1-E2. We tested different
initialization of the Bayesian Optimization, and presented the results of
optimized design parameters and free buckling shape over 100 iterations.

Search trajectory in the SSIM space: The evolution of the design solutions via
iterations to create a target peanut shape is demonstrated in Fig. 4A and
Movie S1.The red dots in Fig. 4A plots the variation of SSIM over iterations.
The maximum SSIMs before certain iteration steps are connected via the blue
line. The maximum SSIM improves significantly from 0.7 (in the first 10 random
searches) to around 0.91 over time. Some red dots are clustered near the
optimal blue line, while some others scatter around the blue line. This
suggests that during the Bayesian optimization, the model explores and
exploits the design space interchangeably. It will not only locally exploit
optimal solutions (such as between iteration 78 and 79 in Fig. 4B, with
materials concentrated to the left and right sides), but also explore unknown
design space (such as iteration 81 and 82 in Fig. 4C. More materials cover the
top and bottom side of the planar structure).

Multiple optimal solutions give SSIM around 0.9: Fig. 4B shows that multiple
latent feature combinations can lead to similar 3D deformed shapes. Even
though these optimal patterns look different, they share the common
characteristics of having more material concentrated near the left and right
ends of the kirigami. The regions with more material compress the bare regions
in between where the stiffness is lower, and such a compression leads to a
ridge in the middle.

Interpolated kirigami patterns are often required to better approximate target
3D shapes: We also explored the generalization capability of VAE and the role
it plays in optimization. For example, the optimal design shown in Fig. 4B2 is
an interpolated pattern generated by the VAE, and not in the patterns used to
train it. We find the nearest neighbours of this optimal generated pattern
using K-nearest neighbors algorithm [36]. If we replace the generated pattern
with its neighbours in in Fig. 4D, we find a decrease in the SSIM. This
suggests the importance of interpolation capability in the latent space to
achieve more accurate design solutions.

The optimal solutions can be different to the initializations: We repeat the
optimization process several times with different initializations (described
in SI Text). For different initializations that lead to similar values of SSIM
(around 0.91) over 100 iterations, the optimized kirigami patterns are
different from each other, as presented in Fig. 4E1 and E2. However, the
optimized prestretch and the radius are close to each other. This suggests
that there exists a region of optimal strain-mismatch and structural size that
leads to the target 3D shape.

Figure 5: Inverse design framework: manufacturing programmable soft 3D
structures. We illustrate the end to end design and manufacturing process.
Target 3D shapes: the first column shows five different 3D shapes (a peanut,
pringle, ship hull, pyramid, flower) input to the optimization algorithm.
Height-coded 2D image representations: the second column shows distribution of
height for the target 3D shape. Optimal design parameters: the third column
presents the optimized kirigami patterns, prestretch, and radius selected
after 100 iterations. Height-coded 2D image representations of the optimal
predicted shapes: the fourth column presents the height distribution of
predicted 3D topology obtained from finite element simulation using the
optimal design parameters. The SSIM are approximately in 0.9. The predicted 3D
shapes: the fifth column shows the predicted 3D shape in simulation. The
maximum heights for the 3D images are shown in the bottom. Experimentally
realized morphing 3D soft structures with desired shapes: the sixth column
shows the experimental result of the 3D topology given the optimal parameters
and following the manufacturing setup illustrated in Fig. 3

We also use the same class of input kirigami patterns to create a pringle and
a ship hull-like 3D shape, presented in Fig. 5B and C. The optimized SSIM over
100 iterations is as high as 0.92 and 0.90, respectively. While the maximum
SSIM over 10 random searches is at 0.70 and 0.73, respectively. This suggests
the importance of strategic search in creating target shapes. The largest
errors are found to be near the sharp edges, where the local curvature has
sharp changes.

###  .2 Target 3D shapes with four-fold radial symmetry

Inspired by flowers in nature and man-made pyramids, we further aim to create
two corresponding deformed soft structures. Both targets are composed of
shapes of four-fold radial symmetry. For example, for each petal in a flower,
the structure is bilaterally symmetric, where each half is a mirror image of
the other half. The generation of such radially symmetric kirigami patterns
are described in SI Text. These kirigami patterns can also be represented by
six latent features, without sacrificing the high reconstruction accuracy
(measured by SSIM $\approx$ 0.99). As shown in Fig. 5, we get very high SSIM.
For the flower like shape, the maximum SSIM increases from 0.81 in the first
10 random combinations to 0.92 within 100 iterations. The maximum SSIM
slightly increases from 0.83 in the first 10 random searches to 0.87 for the
pyramid shape.

###  .3 Experimental validation

We carried out the entire end-to-end design process involving five shapes,
shown in Fig. 5. We picked one optimal design for each target shape from our
algorithm. Next we used these optimal design parameters to manufacture the
corresponding 3D structures. As Fig. 5 shows, the target, the predicted and
the manufactured 3D shapes have the same structural forms, validating our
design process. Given such structural similarity, the maximum height ($H$) of
the 3D structure is an easily measured metric for comparison. As shown in the
fifth and the sixth columns of Fig. 5, the $H$ values are in good agreement
between the predicted and manufactured 3D structures.

One can further explore the physics of the kirigami patterns discovered by our
design framework. For example, for the four-fold radial symmetry involving the
flower and the pyramid target shapes, the optimal kirigami patterns are shown
in Fig. 5E and D. The kirigami pattern for the flower shape has materials
concentrated in the four lobes and has material removed in the center. This
can explain the final shape: the bending of the areas covered with kirigami
creates a ridge of slightly lower height in the middle, giving rise to the
four lobed flower pattern. In contrast, the pyramid shape with the maximum
height in the center requires more material concentrated in the center, as
shown in Fig. 5D.

Figure 6: Scaling law and kirigami invariance: the same kirigami can be used
to design target shapes with varying size parameters $H/R$ by selecting
prestretch and radius of the composite structure A. Distribution of $H/R$ with
the variation of size ratio $R/t$ and the amount of applied prestretch
$\epsilon$. The star symbol denotes the optimal size and prestretch
combination that gives the target pyramid shape. The black dots denote the
regions where the size and prestretch is perturbed around the optimal point.
B. Typical free buckling shapes that correspond to the black dots in A. Both
the height-coded 2D image representation and the corresponding 3D structure
are presented. C. A normalized shape parameter $Ht/R^{2}$ as a function of
normalized radius $R/t$ at different amount of prestretch. The line horizontal
dashed lines indicates the scaling prediction: $H/R$ scales linearly with
$R/t$

## Discussion

### Kirigami invariance: Scaling Law

We address the problem of manufacturing the same shape but with different
scales (e.g., pyramids with different heights), without having to solve for
optimal design parameters for each scale. Recall that designing a target 3D
shape requires around hundreds of optimization steps, each involving an
expensive forward simulation step. Thus, if we can identify a scaling law
involving the radius $R$, prestretch $\lambda$, and the scale of the target
shapes –while fixing the kirigami patterns– then such a scaling law can guide
designers to quickly vary the scale of the deployable 3D soft structures.

For a given material property and kirigami pattern, the maximum normalized
height $H/R$ is affected by prestretch $\lambda=1+\epsilon$, and normalized
radius $R/t$, where $H$ is the maximum height of the free buckling shape, and
$R$ and $t$ represent the radius and thickness of the planar structure,
respectively (see Fig. 6). We derive an analytical relationship by balancing
the stretching-induced energy (pre-buckling) in the bottom layer and the
bending-related energy dominated by the kirigami layer. We find that for a
pyramid like 3D shape, the maximum normalized height $H/R$ scales with
$\frac{R}{t}\sqrt{\frac{E_{s}}{E_{k}}}{\frac{\epsilon}{1+{\epsilon}}}$ (see SI
Text), where $E_{s}$ and $E_{k}$ are the Young’s moduli of the substrate and
the kirigami, respectively. Such a relationship can directly extend the design
of soft composite structures from one maximum height to different heights by
changing the original radius of the composite structures and prestretch.

As an example, we choose the optimized kirigami pattern in Fig. 5D that leads
to a pyramid shape, and investigate how the normalized height is scaled with
the variation of normalized radius and prestretch. Fig. 6A shows the maximum
normalized height $H/R$ as a function of normalized radius $R/t$ and
prestretch. The star symbol in Fig. 6A denotes the optimal combination of
prestretch and radius that leads to the target pyramid shape over 100
iterations. The black dots are the selected examples where the radius and
prestretch are perturbed around the optimal values. If we want to design
taller pyramids, we can slightly increase the normalized radius from Fig. 6B1
to Fig. 6B2. When we slightly increase the prestretch around the optimal point
from Fig. 6B1 to Fig. 6B3, the height of the 3D shape also increases, but the
radius of the deformed geometry decreases. Similar phenomena can be observed
by varying from Fig. 6B2 to Fig. 6B4.

The increase of $H/R$ is proportional to the increase of the normalized radius
$R/t$ in Fig. 6C. This agrees with the scaling analysis of elastic energies.
When the prestretch is larger, the stretching energy in the deformed (i.e.
post-buckling) configuration becomes larger, which causes $H/R$ to deviate
from the linear relationship with $R/t$. As $R/t$ decreases, the more
important the post-buckling stretching energy becomes, but such stretching
energy is not considered in the current analytical derivation, which is based
on balancing only the pre-buckiling stretching energy.This suggests that the
coupling between stretching and bending are important in these conditions,
which needs more complicated models beyond the scaling analysis to explain the
coupling and hence the shape variation. Similar phenomena are also found for
other target shapes, such as flower and peanut, and more detailed discussion
can be found in SI Text.

### Importance of symmetry

In this section, we aim to demonstrate that an arbitrarily constructed large
class of kirigami patterns may still be limited in forming target 3D shapes.
The free buckling shape assumed by the symmetry-constrained shapes and
unidirectional strips, which doesn’t guarantee the reflectional or four-fold
radial symmetry, are compared. Fig. 7A shows that even though the
unidirectional strips are able to predict the ship hull like shape with
reasonable accuracy, it fails to predict the flower shape in Fig. 7B. The
optimized kirigami patterns and the corresponding 3D shapes after 100
iterations are also compared with the target shape in Fig. 7. Even though the
VAE models can generate new patterns, with mixed combination of multi-
directional strips, few of the generated patterns are exactly radially
symmetric. This causes the free buckling shape to bend more in a preferred
direction. The desired reflectional symmetry of a ship hull like target shape
is also not guaranteed using the combinations of unidirectional strips. The
comparisons demonstrate the importance of constraining the design space with
appropriate symmetric property at the very beginning of design process.

Figure 7: Comparison of the optimal kirigami patterns and height distribution
of assumed 3D shape using different optimization methods. Target shapes: The
first column shows the distribution of height for the target shapes A. Ship
hull B. Flower shape. The optimizations are conducted using four approaches.
VAE+Bayesian optimization with symmetric patterns: the second column presents
the optimal solution found in 100 iterations using VAE+Bayesian optimization
with symmetric patterns. VAE+Bayesian optimization with unidirectional strips:
the third column shows the optimal solution in 100 iterations using
VAE+Bayesian optimization with unidirectional strips. Genetic algorithm with
symmetric patterns: the fourth column shows the optimal solution in 200
iterations using genetic algorithm with symmetric patterns. Genetic algorithm
with unidirectional strips: the fifth column shows the optimal solution in 200
iterations using genetic algorithm with unidirectional strips.

### The advantage of proposed VAE and Bayesian Optimization combined approach
over evolutionary algorithms

We also compare the optimization results from the proposed framework with that
using genetic algorithm (GA), which is a standard evolutionary algorithm (see
SI Text). In the evolutionary approach, we set the variables to be 13 binary
integers representing the presence of kirigami materials in the 13 divided
regions in Fig.S1 A2, B2, and C2 . Another discrete integer is used to
indicate the amount of rotation applied to the kirigami patterns, and two
other continuous variables are introduced to represent the radius of the
structure, and the amount of prestretch to be applied. The population size in
the standard algorithm is chosen as 10, and the maximum number of iterations
is set as 20. This means that we need to perform 200 numerical simulations per
design; recall that this is twice the number of forward simulations used in
our design framework, thus giving the GA approach a computational advantage.

Fig. 7 compares the resucolt of applying genetic algorithm to the
unidirectional or symmetric shapes.The comparison suggests that the proposed
VAE and Bayesian Optimization combined approach can search the optimal
combinations of kirigami and prestretch faster, and can achieve 3D deformed
shapes closer to target shapes, compared to the evolutionary search. Without
the VAE-aided reduced dimensional and continuous kirigami search space, the
standard evolutionary optimization is limited to searching only the discrete
design space. The comparisons between the third and the fifth column of Fig.
7A and B demonstrate that the incapability of generating new kirigami patterns
beyond the input dataset restricts the conventional optimization approaches
from interpolating between the candidate kirigami patterns to achieve desired
free buckling shapes.

## Conclusions

We numerically and experimentally explored the capability of deforming a
planar composite structure to target 3D shape via kirigami cutting and strain-
mismatch. The design space is very high dimensional to be optimized directly.
We formulated a VAE and active learning combined approach to tackle the design
challenges. A VAE is used to represent originally high-dimensional design
variables to a much lower dimensional continuous search space. The Bayesian
optimization is then conducted to quickly obtain multiple optimal design
solutions that achieve similar target free buckling shapes, ranging from
shapes inspired by a peanut to a pyramid. We found that the nonlinear
interplay between the strain mismatch, size of the composite structure, and
the kirigami patterns strongly affect the free buckling shapes. We also
studied the effect of imposing symmetry constraints on the machine learning-
aided design results. The comparison of the results with and without symmetry
constraints demonstrates that the symmetry constraints at the beginning of the
machine learning process are important in better approximating the target
shapes. A scaling law is used to guide scaling of the target shapes from one
size to another, without having to search for the optimal design parameters.
We also discussed the advantage of the proposed framework over traditional
approaches, such as genetic algorithm. The proposed framework accelerates the
design of a series of shape morphing, fully soft composite structures from
weeks and months of running millions of simulations to a few hours of
strategically examining around 100 examples.

The inverse design method can provide a systematic way to solve a variety of
form finding problems not limited to soft kirigami structures, but also to the
manufacturing of gridshell [37, 38] and compressive buckling-induced 3D
architectures using micro ribbons [39], which can find applications in a range
of areas including soft robotics, additive manufacturing, and architecture. In
our future work we plan to improve fabrication accuracy even further by
potentially incorporating controls that can induce local deformations into our
planar-only manufacturing platform.

## Materials and methods

### Bayesian optimization

Bayesian optimization is a sample-efficient approach for solving a wide range
of global optimization problems. See Ref. [28] for details; a short summary
follows. This approach aims to solve an optimization, expressed as

$\theta^{*}=\underset{\theta}{\mathrm{argmax}}f(\theta),$ (1)

where $f$ is a black box model which is expensive to evaluate. The function is
approximated by a Gaussian process model. A Gaussian process prior using a
Matern similarity kernel with homoscedastic noise is selected as the
functional prior $\Phi(\theta)$. Given the available data $D$, the posterior
distribution of the parameters is computed via the Bayes’ theorem,

$\Phi(\theta)=\frac{\Phi(D|\theta)\Phi_{prior}(\theta)}{\Phi(D)}$ (2)

Based on the current posterior distribution, the acquisition function $EI$ is
selected as the

$EI(\theta)=\mathbb{E}[u(\theta)]=\mathbb{E}[[f(\theta)-f(\theta_{t}^{*})]^{+}]$
(3)

where $\mathbb{E}$ is computing the expectation, $\theta_{t}^{*}$ is the best
point observed so far. When $f(\theta)>f(\theta_{t}^{*})$,
$u(\theta)=f(\theta)-f(\theta_{t}^{*})$. While when $f(\theta)\leq
f(\theta_{t}^{*})$, $f(\theta_{t}^{*})=0$. The next data to be sampled
$\theta_{t+1}$ is selected such that the acquisition function is maximised,
i.e., $\theta_{t+1}=\underset{\theta}{\mathrm{argmax}EI(\theta)}$. Then, the
Gaussian process and acquisition functions are updated. Such a process
iterates several times, until the optimal solution converges. The algorithm is
implemented using the Scikit-Optimize package [40].

### Implementation details of finite element simulation

In the finite element simulations, the substrate and kirigami layers were
meshed using 3-node triangular shell element. The substrate layer is divided
into 2772 triangle elements. The material properties for the substrate and
kirigami layer are listed in Table LABEL:table:matparam.

Table 1: Material parameters for the substrate and kirigami layer Parameter | Value
---|---
$C_{1}^{\textrm{s}}$ | 2.4 KPa
$C_{2}^{\textrm{s}}$ | 23.4 KPa
Thickness $t^{\textrm{s}}$ | 1.1 mm
$C_{1}^{\textrm{k}}$ | -2.6 KPa
$C_{2}^{\textrm{k}}$ | 185.8 KPa
Thickness $t^{\textrm{k}}$ | 1.4 mm
Outer radius of the circular substrate | 3 cm

##  I Data deposition

We have created a repository at https://github.com/StructuresComp/inverse-
kirigami containing the matlab code for generating finite element simulations
and python code for performing machine learning.

##  II Acknowledgement

We thank Shyan Shokrzadeh and Vishal Kackar for their assistance on
experiments. Following research grants are gratefully acknowledged: NSF
(CMMI-2053971) for L.M., M.M., V.R., and M.K.J.; and NSF (CAREER-2047663,
CMMI-2101751) for M.K.J.

## References

  * [1] JW Boley, et al., Shape-shifting structured lattices via multimaterial 4d printing.  Proceedings of the National Academy of Sciences 116, 20856–20862 (2019).
  * [2] Y Tang, J Yin, Design of multifunctional soft doming actuator for soft machines.  Advanced Materials Technologies 3, 1800069 (2018).
  * [3] MA Dias, et al., Kirigami actuators.  Soft matter 13, 9087–9092 (2017).
  * [4] L Tomholt, et al., Tunable infrared transmission for energy-efficient pneumatic building façades.  Energy and Buildings 226, 110377 (2020).
  * [5] J Panetta, M Konaković-Luković, F Isvoranu, E Bouleau, M Pauly, X-shells: A new class of deployable beam structures.  ACM Transactions on Graphics (TOG) 38, 1–15 (2019).
  * [6] S Abdelmohsen, et al., A multi-physics approach for modeling hygroscopic behavior in wood low-tech architectural adaptive systems.  Computer-Aided Design 106, 43–53 (2019).
  * [7] S Xu, et al., Assembly of micro/nanomaterials into complex, three-dimensional architectures by compressive buckling. Science 347, 154–159 (2015).
  * [8] J Kim, JA Hanna, M Byun, CD Santangelo, RC Hayward, Designing responsive buckled surfaces by halftone gel lithography.  science 335, 1201–1205 (2012).
  * [9] AE Forte, et al., Inverse design of inflatable soft membranes through machine learning.  Advanced Functional Materials 32, 2111610 (2022).
  * [10] S Felton, M Tolley, E Demaine, D Rus, R Wood, A method for building self-folding machines.  Science 345, 644–646 (2014).
  * [11] M Pezzulla, GP Smith, P Nardinocchi, DP Holmes, Geometry and mechanics of thin growing bilayers.  Soft matter 12, 4435–4442 (2016).
  * [12] C Baek, AO Sageman-Furnas, MK Jawed, PM Reis, Form finding in elastic gridshells.  Proceedings of the National Academy of Sciences 115, 75–80 (2018).
  * [13] M Liu, L Domino, D Vella, Tapered elasticæ as a route for axisymmetric morphing structures.  Soft Matter 16, 7739–7750 (2020).
  * [14] Z Fan, et al., Inverse design strategies for 3d surfaces formed by mechanically guided assembly.  Advanced Materials 32, 1908424 (2020).
  * [15] WM Van Rees, E Vouga, L Mahadevan, Growth patterns for shape-shifting elastic bilayers.  Proceedings of the National Academy of Sciences 114, 11597–11602 (2017).
  * [16] X Dang, et al., Inverse design of deployable origami structures that approximate a general surface.  International Journal of Solids and Structures 234, 111224 (2022).
  * [17] J Wang, J Suo, A Chortos, Design of fully controllable and continuous programmable surface based on machine learning.  IEEE Robotics and Automation Letters 7, 549–556 (2021).
  * [18] M Moshe, et al., Kirigami mechanics as stress relief by elastic charges.  Physical review letters 122, 048001 (2019).
  * [19] Y Zhang, et al., A mechanically driven form of kirigami as a route to 3d mesostructures in micro/nanomembranes.  Proceedings of the National Academy of Sciences 112, 11757–11764 (2015).
  * [20] DP Holmes, Elasticity and stability of shape-shifting structures.  Current opinion in colloid & interface science 40, 118–137 (2019).
  * [21] L Stein-Montalvo, P Costa, M Pezzulla, DP Holmes, Buckling of geometrically confined shells.  Soft Matter 15, 1215–1222 (2019).
  * [22] J Zavodnik, Y Wang, W Yan, M Brojan, MK Jawed, Soft kirigami composites for form-finding of fully flexible deployables.  arXiv preprint arXiv:2301.06597 1, 1 (2023).
  * [23] Z Liu, D Zhu, SP Rodrigues, KT Lee, W Cai, Generative model for the inverse design of metasurfaces.  Nano letters 18, 6570–6576 (2018).
  * [24] PZ Hanakata, ED Cubuk, DK Campbell, HS Park, Forward and inverse design of kirigami via supervised autoencoder.  Physical Review Research 2, 042006 (2020).
  * [25] R Danhaive, CT Mueller, Design subspace learning: Structural design space exploration using performance-conditioned generative modeling.  Automation in Construction 127, 103664 (2021).
  * [26] Y Mao, Q He, X Zhao, Designing complex architectured materials with generative adversarial networks.  Science advances 6, eaaz4169 (2020).
  * [27] PZ Hanakata, ED Cubuk, DK Campbell, HS Park, Accelerated search and design of stretchable graphene kirigami using machine learning.  Physical review letters 121, 255304 (2018).
  * [28] PI Frazier, A tutorial on bayesian optimization.  arXiv preprint arXiv:1807.02811 - (2018).
  * [29] JD Paulsen, et al., Curvature-induced stiffness and the spatial variation of wavelength in wrinkled sheets.  Proceedings of the National Academy of Sciences 113, 1144–1149 (2016).
  * [30] TC Shyu, et al., A kirigami approach to engineering elasticity in nanocomposites through patterned defects.  Nature materials 14, 785–789 (2015).
  * [31] A Rafsanjani, K Bertoldi, Buckling-induced kirigami.  Physical review letters 118, 084301 (2017).
  * [32] B Burger, et al., A mobile robotic chemist.  Nature 583, 237–241 (2020).
  * [33] E Ahmed, et al., A survey on deep learning advances on different 3D data representations.  arXiv preprint arXiv:1808.01462 -, – (2018).
  * [34] H Zhao, O Gallo, I Frosio, J Kautz, Loss functions for image restoration with neural networks.  IEEE Transactions on computational imaging 3, 47–57 (2016).
  * [35] M Smith, ABAQUS/Standard User’s Manual, Version 6.9.  (Dassault Systèmes Simulia Corp, United States), (2009).
  * [36] F Pedregosa, et al., Scikit-learn: Machine learning in python.  the Journal of machine Learning research 12, 2825–2830 (2011).
  * [37] L Qin, W Huang, Y Du, L Zheng, MK Jawed, Genetic algorithm-based inverse design of elastic gridshells.  Structural and Multidisciplinary Optimization 62, 2691–2707 (2020).
  * [38] F Isvoranu, J Panett, T Chen, E Bouleau, M Pauly, X-shell pavilion: A deployable elastic rod structure in Proceedings of IASS Annual Symposia.  (International Association for Shell and Spatial Structures (IASS)), Vol. 2019, pp. 1–8 (2019).
  * [39] H Luan, et al., Complex 3d microfluidic architectures formed by mechanically guided compressive buckling.  Science Advances 7, eabj3686 (2021).
  * [40] T Head, M Kumar, H Nahrstaedt, G Louppe, I Shcherbatyi, Scikit-optimize/scikit-optimize.  (version 0.8. 1) - (2020).

Supplementary Information

### Creating symmetric kirigami patterns

For example, when the 3D structure has two symmetry axes, we divide the 2D
structure into four quadrants. According to the symmetric property of target
3D structure, the ”planar design space” reduces to the regions in the first
quadrants. In the first quadrant, the kirigami patterns are created by
randomly combining the 13 independent pieces, as shown in Figure 8 A1, B1 and
C1. Then, in case of reflectional symmetry, the kirigami patterns are created
by reflecting the image in the first quadrant with respect to the two symmetry
axes, as shown in Figure 8 B3.

While in case of structures that are of four-fold radial symmetry, the planar
design space further reduces to the half of the first quadrant. The other half
can be directly created by the mirror image in Figure 8 C3. The final kirigami
cuts are created by rotating the patterns of the first quadrants by 90, 180,
and 270 degrees, as presented in Figure 8 C4.

### Augmenting the kirigami design space via rotation

Rotation is a common way to introduce more diverse patterns for machine
learning training. Once we create a certain set of kirigami cuts in the design
space, we further rotate the kirigami cuts in the clockwise and
counterclockwise directions two to four times, and only keep the rotated
patterns inside the design space, as shown in Figure 8 A2, B2, and C2. For the
class of kirigami with unidirectional strips, we rotate the strips four times,
which are 72, 144, 216, 288 degrees in the clockwise direction. For the class
of the kirigami patterns of reflectional symmetry, the kirigami cuts are
rotated 7.5 and 15 degrees in the counter-clockwise direction, and are rotated
in 4.5 and 12 degrees in the clockwise direction. Finally, for the class of
the patterns with four-fold radial symmetry, the cuts in the half of the first
quadrants are rotated in 7.5 degrees in the counter-clockwise direction, and
in 4.5 degrees in the clockwise direction.

Figure 8: Key steps for creating kirigami patterns for A. Unidirectional
strips. B. Kirigami patterns that have reflectional symmetry. C. Kirigami
patterns that have four-fold radial symmetry

### Variational Autoencoder

We assume that each observed data point $y_{i}$ is generated in a nonlinear
fashion with latent variable $z_{i}$, whose joint probability density can be
expressed as,

$p(y,z)=\prod_{i=1}^{N}p_{g}(y_{i}|z_{i})p(z_{i})$ (4)

The full variational lower bound to be optimized include a reconstruction loss
(Binary Cross Entropy Loss) and a Kullback-Leibler divergence term, which is
written as,

$ELBO=\mathbb{E}_{z\sim{q(z|y)}}{\log p(y|z)}-D_{KL}(q(z|y)||p(z))$ (5)

The approximate posterior given by the inference network for image $y$ is
$q(z|y)=N(\mu_{z}(y),\sigma_{z}^{2}(y))$, which follows a normal distribution
with mean $\mu_{z}(y)$ and standard deviation $\sigma_{z}(y)$. The prior
distribution is given as a normal distribution $p(z)=N(0,I)$.

### Structure of the Variational AutoEncoder

The weights in the encoder and decoder are optimized simultaneously to
minimize the ELBO loss function, including the KL divergence term and the
reconstruction loss. The encoder and decoder each has two hidden layers with
128 neurons per layer. The activation function for each neuron is chosen as
the $tanh()$ function. While the activation function for the output is chosen
as $sigmoid()$ function. The loss are minimized using the Adam optimization,
an extension of Stochastic Gradient Descent with the learning rate set as
0.0005. The input and output for the VAE model are the same binary image,
whose size is 64 by 64.

### Examples of generated kirigami patterns via the Variational AutoEncoder

Figure 9 presents some typical generated images when images of reflection
symmetry (data size of training set is 40955) are input to the Variational
Autoencoder model. These images are created by inputting 10 cluster centers in
a K-Means clustering process into the trained decoder neural network. It can
be found the typical images are also reflection symmetry.

Figure 9: Typical kirigami patterns, reconstructed from the 5 cluster centers
of the latent features using K-Means clustering A. Reflectional symmetry
patterns B. Four-fold radial symmetry patterns

### The effect of initialization on the results for Bayesian optimization

For the peanut shape, we redo the Bayesian optimiaztion with different initial
choices of 10 random searches. We find that the initialization affects the
maximum SSIM over 100 iterations in Figure 10. However, for the
initializations that give similar maximum SSIM (Initialization No. 1, 2, 4),
the amount of optimized prestretch and kirigami radius are found to be
similar, in Figure 10B and C. This suggests that there exists a sweet-spot of
strain-mismatch and structural size that leads to the target 3D shape.

Figure 10: For different initalizations, the Bayesian Optimization gives the
(A) Variation of maximumm SSIM over iterations. (B) Variation of kirigami
radius over iterations. (C) Variation of prestretch over iterations (Target
shape: peanut shape)

### Detailed description of the genetic algorithm

Table 2: Parameters for genetic algorithm Parameter | Value
---|---
Maximum number of iterations | 10
Population size | 20
Mutation probability | 0.1
Elite ratio | 0.01
Crossover type | uniform
Crossover probability | 0.5
Parents portion | 0.3
Crossover type | uniform

The genetic algorithm is a search heuristic relying on biologically inspired
operators such as mutation, crossover and selection. The idea of selection
phase is to select the fittest individuals and let them pass their genes to
the next generation.The “parents” are selected, based on their fitness scores
(e.g. the error between prediction and target). Individuals with high fitness
have more chance to be selected as parents for reproduction. The offsprings
are generated by exchanging the genes of the parents at certain
probability.These generated offsprings are added into the population. To
maintain the diverstiy within the population, a portion of the genes are
replaced by a random value. This process is also called as mutation. A
pracitcal variant of the general genetic algorithm process also involves
allowing the best individuals, determined by the elite ratio,to be directly
carry over to the next generation, unaltered. This helps guarantees that the
solution quality obtained by the genetic algorithm won’t decrease from
generation to the next. The algorithm is implemented using
https://github.com/rmsolgi/geneticalgorithm.

When the genetic algorithm is applied to the kirigami composite designs, the
inputs (i.e. genes) are variables mixed with integers and continous variables.
The fitness function to be minimized in the genetic algorithm is the 1-SSIM,
where SSIM measures the structural similarity between the prediction and the
target. The parameters for the genetic algorithm is listed in Table
LABEL:table:genparam.

### Trade-off between size of the composite and the amount of pre-stretch

If we keep the kirigami fixed, the target shapes of other 3D geometries also
have nonlinear relationships with prestretch and kirigami size, especially
when prestretches are large, as shown in 12 and 13. In general, we can see
that the maximum height $H/R$ increases together with the prestretch and
structure size. The radius of the deformed shape increases with the structure
size, but decreases with larger prestretch. Even though the linear
relationship between $H/R$ and $R/h$ holds well for small prestretch
conditions, as prestretch increases, the 3D shape can be very different from
the deformation under small prestretch. The actual relation between $H/R$ and
$R/h$ starts deviating from the linear relationship. For instance, for the
kirigami pattern that creates the targeted peanut shape, if we decrease the
prestretch to a small enough value, the 3D shape also turns into a pringle
shape, rather than maintaining the peanut shape (13).

### Simulation of prestretch release via temperature reduction

To simulate the prestretch release, we apply an isotropic expansion.

$\displaystyle\epsilon_{sim}=\alpha(\theta-\theta^{0})=\alpha\delta\theta$ (6)

where $\alpha$ is the thermal expansion coefficient, which is set to 1,
$\epsilon_{sim}$ is the strain induced by the temperature variation. $\theta$
and $\theta_{0}$ are the current and initial temperature, respectively. The
temperature variation $\delta\theta$ is relatd to the prestretch in the
experiments $\epsilon$ via

$\displaystyle\epsilon=\frac{\delta\theta}{1-\delta\theta}$ (7)

### Analytical analysis of the effect of design parameters on the free
buckling shape

For small strain, we calculate the stretching-induced energy for the substrate
layer as,

$\displaystyle{\\!\\!\text{E}\,\,}_{s}=\frac{1}{2}{\pi{R^{2}}}C_{s}\delta\theta^{2}$
(8)

where $C_{s}=E_{s}t/(1-\nu^{2})$.$t$ is the thickness of the plate, which is
around $1mm$ in the experiment.

We assume that after the deformation, the structure is bending dominated, and
the curvature is about the same everywhere. The bending energy can be
calculated as,

$\displaystyle{\\!\\!\text{E}\,\,}_{b}$
$\displaystyle=\frac{1}{2}{\pi{R^{2}}}D_{s}\kappa^{2}+\frac{1}{2}n{\pi{R^{2}}}D_{k}\kappa^{2}\approx\frac{1}{2}n\pi{R^{2}}D_{k}\kappa^{2}$
(9)

where $n$ is the fraction of the area covered by the kirigami layer, which is
a factor between 0 and 1.
$D_{s}=\frac{{E_{s}{t^{3}}}}{{12\left({1-{\nu^{2}}}\right)}}$ and
$D_{k}=\frac{{E_{k}{t^{3}}}}{{12\left({1-{\nu^{2}}}\right)}}$ are the bending
rigidity of the substrate and kirigami layer, respectively.

We further assume that the deformed shape is close to a spherical cap shape
with maximum height $H$, and has approximately constant curvature everywhere,
as shown in Figure 11. The length of the curve $R$ can be approximated as,

Figure 11: Sketch of a spherical cap with height $H$

$\displaystyle R^{2}=\rho^{2}-(\rho-H)^{2}+H^{2}$ (10)

where $\rho$ is the radius of the hemisphere, which is related to the
curvature $\kappa$ via $\kappa=1/\rho$.

Hence, the maximum height of the shape is related to the curvature via,

$\displaystyle H=\frac{R^{2}}{2}\kappa$ (11)

Hence, equalizing the bending energy and stretching energy, and combining
equation (5) gives,

$\displaystyle\frac{H}{R}=\sqrt{\frac{3E_{s}}{nE_{k}}}\delta\theta\frac{R}{t}\approx\sqrt{\frac{3E_{s}}{nE_{k}}}{\frac{\epsilon}{1+{\epsilon}}}\frac{R}{t}$
(12)

$\displaystyle\frac{H}{R}\approx\sqrt{\frac{3E_{s}}{nE_{k}}}{\frac{\epsilon}{1+{\epsilon}}}\frac{R}{t}$
(13)

Figure 12: (A) Distribution of $H/R$ with the variation of size ratio $R/h$
and the amount of applied prestretch $\epsilon$. The star symbol denotes the
optimal size and prestretch combination that gives the target pyramid shape.
The black dots denote the regions where the size and prestretch is perturbed
around the optimal point. (B) Typical free buckling shapes that correspond to
the black dots in (A). (C) Normalized $Hh/R^{2}$ as a function of normalized
radius $R/h$ at different amount of prestretch. The line horizontal dashed
lines indicates the scaling prediction. (Target shape: flower shape) Figure
13: (A) Distribution of $H/R$ with the variation of size ratio $R/h$ and the
amount of applied prestretch $\epsilon$. The star symbol denotes the optimal
size and prestretch combination that gives the target pyramid shape. The black
dots denote the regions where the size and prestretch is perturbed around the
optimal point. (B) Typical free buckling shapes that correspond to the black
dots in (A). (C) Normalized $Hh/R^{2}$ as a function of normalized radius
$R/h$ at different amount of prestretch. The line horizontal dashed lines
indicates the scaling prediction. (Target shape: peanut shape)