# Evaluation of Synthetic Datasets for Conversational Recommender Systems

Harsh Lara
Google Research
<EMAIL_ADDRESS>
Manoj Tiwari
Google Research
<EMAIL_ADDRESS>

###### Abstract

For researchers leveraging Large-Language Models (LLMs) in the generation of
training datasets, especially for conversational recommender systems - the
absence of robust evaluation frameworks has been a long-standing problem Peng
et al. (2017). The efficiency brought about by LLMs in the data generation
phase is impeded during the process of evaluation of the generated data, since
it generally requires human-raters to ensure that the data generated is of
high quality and has sufficient diversity. Since the quality of training data
is critical for downstream applications, it is important to develop metrics
that evaluate the quality holistically and identify biases. In this paper, we
present a framework that takes a multi-faceted approach towards evaluating
datasets produced by generative models and discuss the advantages and
limitations of various evaluation methods.

##  1 Introduction

Conversational recommender systems consist of several machine learning models
working together under the hood to perform complex tasks while providing users
with a seamless and uncomplicated natural language interface. Even the most
primitive conversational recommenders consist of at least four different
components: a NLU unit, a NLG unit, a dialogue policy unit and a retrieval and
ranking system. Recently, the performance of all of these components have
shown promising results using deep neural networks, which require an abundance
of training data to train on.

The advent of pre-trained language models (BERT, GPT-3, etc.) have made it
possible for natural-language applications to achieve admirable performance
using zero-shot or few-shot learning. However, the bulk of complex tasks (e.g.
intent prediction, summarization, providing item recommendations etc.) still
require that these pre-trained models be finetuned for meeting more composite
objectives.

Data collection for fine-tuning tasks often turns out to be the bottleneck
that is the differentiator between successful and unsuccessful machine
learning applications. Data generation using LLMs (e.g. GPT-3, LaMDA, PALM)
Brown et al. (2020); Thoppilan et al. (2022); Chowdhery et al. (2022) is a
promising approach towards treating the data scarcity problem. If we can
leverage LLMs using priming along with pre-programmed dialogue policies to
generate synthetic data that can replicate human interactions with
conversational agents, it can enable a plethora of fine-tuned applications.

Evaluating the quality of synthetic data generated with the help of LLMs is a
challenging task. With the help of LLMs, it is possible to generate thousands
of examples of interactions between the conversational recommender and a
simulated user, however - it is difficult to ensure that the synthesized data
is of good quality. Researchers generally hold synthetically generated data at
a higher bar for quality as compared to human-generated data. The reason for
this is simple - it is easy for biases to creep into data generation when
using LLMs instead of humans. Additionally, it is non-trivial to ensure that
synthesized data is diverse and representative of the gamut of human
interactions with a conversational agent.

Figure 1: Example of a conversation between a user and a conversational movie-
recommender system

##  2 Conversational Recommender Systems

Conversational recommender systems (or conversational recommenders) are
conversational agents which are tasked with engaging in task-oriented, multi-
turn dialogue with users in order to provide useful recommendations to them
Jannach et al. (2020). During such dialogue, the system would elicit user
preferences, provide explanations for recommended items, and process user’s
feedback on the recommendations.

Conversational recommenders are expected to perform better than non-
conversational recommender systems, because of their interactive nature and
their ability to enable high-bandwidth communication with the user by means of
natural language dialogue Rich (1979).

Formally, a conversational recommender can be defined as a conversational
agent (_C_) that produces a combination of natural language response (_S i_)
coupled with a slate of recommended items (_R i_) to the user, after the _i
th_ user utterance (_U i_). A conversation session consists of a sequence of
alternating user and system turns, along with recommendations accompanying
each system turn as shown in Figure 1.

##  3 Datasets for Conversational Recommenders

The objective of conversational recommender systems is two-fold:

  1. 1.

Efficiently recommend items to the user that they would find useful and
interesting, given their preferences at that time.

  2. 2.

Create exploration opportunities for the user, so that they are able to
discover items that they wouldn’t have otherwise found without the assistance
of a conversational recommender system.

Therefore, while developing fine-tuning datasets, it is critical that the data
exhibits good quality with respect to both of these goals. If we were to
optimize for only efficiency of the recommender, we would risk ending up with
a dataset that is degenerate - for example, a dataset in which the recommender
system directly suggests the top-ranked item which matches the user’s
interests in the first turn of the conversation, without creating any
exploration opportunities for the user. If we were to optimize only for
exploration opportunities for the user, we would end up with a dataset
consisting of long-winding conversations that don’t necessarily exhibit
satisfactory efficiency by the conversational recommender.

##  4 Evaluating Conversational Recommenders Datasets

Due to the reasons discussed previously, datasets generated for fine-tuning
conversational recommenders should be evaluated by using a multi-dimensional
approach.

We propose two main criteria for the evaluation of dataset quality:

  1. 1.

_Realism_ : The ability of the dataset to imitate the conversational ability
of real humans.

  2. 2.

_Diversity_ : The ability of fine tuning datasets to exhibit the gamut of
interactions real humans would have had with the recommender system.

###  4.1 Realism

When simulating user interactions with the system, having the conversation
sessions closely mimic real human interactions with the agent is a desirable
property. Intuitively, the more closely we model real human conversations in
our training data generation, the more likely we are to perform better at
downstream tasks such as intent prediction and item reranking.

The natural language utterances from the user enable the conversational
recommender to build a useful model of the user’s preferences and provide
valuable recommendations. Similarly, the natural language output from the
conversational recommender enable the user to have a better understanding of
the system’s suggestions and stimulate a rich back and forth between both
entities. Therefore, during the data generation phase, it is imperative to
produce a high quality of natural language dialogue during the sessions.

Ultimately the goal of the conversational recommender is to model the user’s
preferences better and provide useful recommendations to users while enabling
a high bandwidth interaction between the user and the system. With high
quality dialogue, we expect recommendations produced at the later turns of a
session to be much better than the recommendations produced during the initial
turns of the conversation.

####  4.1.1 Using human raters

Realism can be measured by having crowd-workers evaluate a set of randomly
sampled conversations from a data generation run.

For example, if we have a version of a conversation simulator (say, Simulator
K), we would use that to generate _N_ (say, N=50) conversations. All of these
50 conversations can be sent to crowd-workers with the following instructions:

“Here are 50 user conversations with a conversational recommender system”.

Some conversations in this set are human-generated, whereas others are
simulator-generated. For each conversation, identify whether it is human-
generated or simulator-generated.”

Therefore, we would have the following ratio (let’s call it Realism Score)
computed for the evaluation set that is shared with crowd-workers.

$RealismScore_{human}=\frac{h_{human}}{N}$

Where _h_ is the number of conversations identified as human-generated and _N_
is the total number of conversations shared with human raters.

####  4.1.2 Using a machine learned classifier

Another way of evaluating realism is by training a discriminator to classify
conversations into one of two categories: human-generated or simulated.

A low-cost means of training such a discriminator, is to use a human-rated
batch of conversations as the ground-truth. As discussed in the previous
section, we would share a batch simulated conversations and have crowd-workers
provide binary labels for each of them - _simulator-generated_ or _human-
generated_. the labeled dataset produced in this manner would be used as the
training data for the discriminator.

We can run the discriminator on a batch of simulated conversations to obtain a
Realism Score derived along the same idea we have used with human-raters.

$RealismScore_{inferred}=\frac{h_{predicted}}{N}$

Where _h predicted_ is the number of conversations predicted to be human-
generated and _N_ is the total number of conversations used during inference.

###  4.2 Diversity

Simulator Version | HSentiment | HTopic | Entropy Score
---|---|---|---
A | 0.76 | 0.58 | 0.67
B | 0.82 | 0.67 | 0.745
C | 0.66 | 0.93 | 0.795
Table 1: Examples of entropy scores computed for 3 versions of a conversation
simulator.

In the context of synthetic data generation, we define diversity as the
measure of how balanced the distribution of data points (conversations) is
with respect to the vocabulary over our chosen dimensions (e.g. sentiments,
topics, intents, etc.). As noted by Hashimoto et al. (2019) human evaluation,
which is often viewed as the gold standard evaluation, captures quality
perfectly but fails to capture diversity

We propose the measurement of diversity by evaluating the distribution of
data-points (in our case, dialogue sessions) across different pre-defined
dimensions. For each dimension, we would need to specify a vocabulary (classes
that any data-point can be classified into) and a corresponding classifier (a
system that can classify each data-point into one of the previously defined
classes).

For example, a framework to evaluate conversational datasets could include the
following classifiers:

  1. 1.

Sentiment Classifier

    1. (a)

Vocabulary: Negative, Positive.

    2. (b)

Goal: Classify each session into one of the predefined sentiment categories,
based on the perceived user sentiment.

  2. 2.

Topic Classifier

    1. (a)

Vocabulary: Automobiles, Cooking, Travel, Technology, Fashion, History.

    2. (b)

Goal: Classify each simulated session into a topic class, based on the main
focus of the discussion.

For any data generation job, we can produce an average entropy score, by using
the individual entropies produced by these classifiers.

For each dimension, we would estimate the entropy of the source (data
generator) by using the standard formula for Shannon entropy Shannon (2001):

$H(X)=\sum_{i=1}^{n}p(x_{i})log_{b}\frac{1}{p(x_{i})}$

We would compute a final _Entropy Score_ by taking a weighted average of all
the individual entropy values.

$EntropyScore=\alpha_{1}H_{1}+\alpha_{2}H_{2}+…+\alpha_{n}H_{n}$

Where Hi is the entropy computed for classifier _i_ and $\alpha$i is the
weight assigned for the corresponding entropy. And $\sum_{i=1}^{n}\alpha$i =
1.

We can quantitatively compare the average entropy scores after running data
generation jobs using different versions of the data generator.

For example, consider the situation where we have multiple versions of a
conversation simulator - versions _A_ , _B_ and _C_ with entropy values listed
in Table 1. We would compute the weighted average of entropy values using
equal weights for both entropy values. In this case, simulator version _C_
will rank the highest on the diversity axis, because of its entropy score.

##  5 Conclusion and Future Work

In this paper, we introduce the challenges associated with generating training
data for conversational recommender systems and proposes a multi-dimensional
approach towards the evaluation of the quality of such datasets synthesized
using LLMs. Since we are only beginning to leverage the capabilities of LLMs,
we expect future work in this domain to introduce more elaborate evaluation
strategies that would enable researchers to apply rigorous quality checks on
datasets synthesized using generative models.

## Limitations

The approaches discussed in this paper have several limitations that could be
addressed by future work on evaluation frameworks. The evaluation framework
proposed in this paper is centered around two dimensions that we have found
useful towards setting up evaluation criteria for finetuning datasets for
conversational recommenders, namely, realism and diversity. However, the
proposed evaluation does not take into account the quality of recommendations,
the performance of the agent in terms of identifying user-preferences,
eliciting high-signal information from the user, etc. The work in this paper
is focused mainly on evaluating the quality of natural language dialogue in
the conversational recommender.

## Ethics Statement

This paper aims to initiate a discussion around responsible evaluation of
datasets generated by harnessing the capabilities of generative language
models. In this work, we bring forth the idea that it is important to
carefully evaluate the quality of training datasets generated using LLMs so
that we can ensure good performance on downstream applications that rely on
that training data.

## Acknowledgements

We would like to acknowledge the help and support of all the members on our
team for reviewing the paper and providing valuable feedback.

## References

  * Brown et al. (2020) Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, Christopher Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei. 2020.  Language models are few-shot learners.  _CoRR_ , abs/2005.14165.
  * Chowdhery et al. (2022) Aakanksha Chowdhery, Sharan Narang, Jacob Devlin, Maarten Bosma, Gaurav Mishra, Adam Roberts, Paul Barham, Hyung Won Chung, Charles Sutton, Sebastian Gehrmann, Parker Schuh, Kensen Shi, Sasha Tsvyashchenko, Joshua Maynez, Abhishek Rao, Parker Barnes, Yi Tay, Noam Shazeer, Vinodkumar Prabhakaran, Emily Reif, Nan Du, Ben Hutchinson, Reiner Pope, James Bradbury, Jacob Austin, Michael Isard, Guy Gur-Ari, Pengcheng Yin, Toju Duke, Anselm Levskaya, Sanjay Ghemawat, Sunipa Dev, Henryk Michalewski, Xavier Garcia, Vedant Misra, Kevin Robinson, Liam Fedus, Denny Zhou, Daphne Ippolito, David Luan, Hyeontaek Lim, Barret Zoph, Alexander Spiridonov, Ryan Sepassi, David Dohan, Shivani Agrawal, Mark Omernick, Andrew M. Dai, Thanumalayan Sankaranarayana Pillai, Marie Pellat, Aitor Lewkowycz, Erica Moreira, Rewon Child, Oleksandr Polozov, Katherine Lee, Zongwei Zhou, Xuezhi Wang, Brennan Saeta, Mark Diaz, Orhan Firat, Michele Catasta, Jason Wei, Kathy Meier-Hellstern, Douglas Eck, Jeff Dean, Slav Petrov, and Noah Fiedel. 2022\.  Palm: Scaling language modeling with pathways.
  * Clark et al. (2021) Elizabeth Clark, Tal August, Sofia Serrano, Nikita Haduong, Suchin Gururangan, and Noah A. Smith. 2021.  All that’s ’human’ is not gold: Evaluating human evaluation of generated text.
  * Gehrmann et al. (2022) Sebastian Gehrmann, Elizabeth Clark, and Thibault Sellam. 2022.  Repairing the cracked foundation: A survey of obstacles in evaluation practices for generated text.
  * Hashimoto et al. (2019) Tatsunori B. Hashimoto, Hugh Zhang, and Percy Liang. 2019.  Unifying human and statistical evaluation for natural language generation.  _CoRR_ , abs/1904.02792.
  * Jannach et al. (2020) Dietmar Jannach, Ahtsham Manzoor, Wanling Cai, and Li Chen. 2020.  A survey on conversational recommender systems.  _CoRR_ , abs/2004.00646.
  * Kumar et al. (2020) Varun Kumar, Ashutosh Choudhary, and Eunah Cho. 2020.  Data augmentation using pre-trained transformer models.  _CoRR_ , abs/2003.02245.
  * Peng et al. (2017) Baolin Peng, Xiujun Li, Lihong Li, Jianfeng Gao, Asli Celikyilmaz, Sungjin Lee, and Kam-Fai Wong. 2017.  Composite task-completion dialogue policy learning via hierarchical deep reinforcement learning.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_. Association for Computational Linguistics.
  * Rich (1979) Elaine Rich. 1979.  User modeling via stereotypes.  _Cognitive Science_ , 3(4):329–354.
  * Shannon (2001) C. E. Shannon. 2001.  A mathematical theory of communication.  _SIGMOBILE Mob. Comput. Commun. Rev._ , 5(1):3–55.
  * Thoppilan et al. (2022) Romal Thoppilan, Daniel De Freitas, Jamie Hall, Noam Shazeer, Apoorv Kulshreshtha, Heng-Tze Cheng, Alicia Jin, Taylor Bos, Leslie Baker, Yu Du, YaGuang Li, Hongrae Lee, Huaixiu Steven Zheng, Amin Ghafouri, Marcelo Menegali, Yanping Huang, Maxim Krikun, Dmitry Lepikhin, James Qin, Dehao Chen, Yuanzhong Xu, Zhifeng Chen, Adam Roberts, Maarten Bosma, Yanqi Zhou, Chung-Ching Chang, Igor Krivokon, Will Rusch, Marc Pickett, Kathleen S. Meier-Hellstern, Meredith Ringel Morris, Tulsee Doshi, Renelito Delos Santos, Toju Duke, Johnny Soraker, Ben Zevenbergen, Vinodkumar Prabhakaran, Mark Diaz, Ben Hutchinson, Kristen Olson, Alejandra Molina, Erin Hoffman-John, Josh Lee, Lora Aroyo, Ravi Rajakumar, Alena Butryna, Matthew Lamm, Viktoriya Kuzmina, Joe Fenton, Aaron Cohen, Rachel Bernstein, Ray Kurzweil, Blaise Aguera-Arcas, Claire Cui, Marian Croak, Ed H. Chi, and Quoc Le. 2022.  Lamda: Language models for dialog applications.  _CoRR_ , abs/2201.08239.


# A Hessian-based energy-optimization study of morphoelasticity

Min Wu<EMAIL_ADDRESS>Department of Mathematical Sciences, Worcester
Polytechnic Institute  Center for Computational Biology, Flatiron Institute

###### Abstract

Morphoelasticity describes a family of nonlinear PDE system, and has been
widely applied to understand growth and morphogenesis in living and nonliving
structures. Most previous studies depend on locally linearizing either the PDE
system or its discretized numerical counterpart. This involves conducting a
linear stability analysis for a trivial solution or utilizing Newton’s
iteration on the force residual equations derived from finite element
approximation, respectively. The paper presents a numerical approach that
frames the search for solutions as the optimization of a discretized energy
functional, similar to Tallinen et al., PRL, 2013. In contrast to Tallinen et
al., we utilize the Hessian matrix of the energy to develop a gradient-descent
algorithm that ensures energy stability, which allows us to search solutions
robustly in different types of bifurcations. The solutions are further refined
locally by Newton’s iterations. We showcase the effectiveness of our approach
by studying shape transitions induced by the growth in the interior layer of a
bilayer annulus, with the outer layer constrained at its outermost boundary.
By investigating three cases with varying elastic-moduli ratios between the
two layers, we illustrate that our new method can not only identify and
characterize supercritical bifurcations leading to smooth wrinkles but also
capture subcritical transitions from smooth solutions to non-smooth crease-
forming solutions. Moreover, it reveals aperiodic solutions, both smooth and
non-smooth, that have not been previously reported.

##  1 Introduction

The growth-elasticity theory (aka morphoelasticity) has been widely applied to
understand folding patterns and muscle-contraction related functions in living
organisms by various research groups (see the first work [19], the recent
review [1] and book [8]). In this theory, the dynamics of volumetric growth is
considered to be much slower than the elastic responses. Thus, we can think of
growth as a proxy to load internal stresses into a continuum domain which is
subject to mechanical equilibrium.

In this theory, one follows the Kröner-Lee approach [13, 14] to decompose the
deformation gradient
$\mathbf{F}={\partial\mathbf{x}}/{\partial\mathbf{X}}:=\nabla_{0}\varphi$ by

$\mathbf{F}=\mathbf{F}_{e}\mathbf{G}$ (1)

where $\mathbf{x}=:\varphi(\mathbf{X},t)$ is the observable motion, mapping a
particle at $\mathbf{X}$ initially to its current position $\mathbf{x}$,
$\mathbf{G}(\mathbf{X},t)$ is the accumulated growth (stretch) tensor, and
$\mathbf{F}_{e}$ is the elastic deformation tensor associated with stress.

When the trajectory of $\mathbf{G}(\mathbf{X},t)$ is prescribed, one can solve
$\varphi(\mathbf{X},t)$ and $\mathbf{F}_{e}=\mathbf{F}\mathbf{G}^{-1}$ through
mechanical equilibrium

$\displaystyle\begin{cases}\nabla_{0}\cdot\bm{\Pi}&=\mathbf{0},\hskip
28.45274pt\bm{\Pi}=J_{g}\frac{\partial
W}{\partial\mathbf{F}_{e}}\mathbf{G}^{-T}\text{ in }\Omega_{0}\\\
\bm{\Pi}\mathbf{N}&=\mathbf{T}\text{ on }\partial\Omega_{0}^{N}\end{cases}$
(2)

which is the Euler–Lagrange equation of the energy functional

$\displaystyle
E=\int_{\Omega_{0}}J_{g}W(\mathbf{F}_{e})d\mathbf{X}-\int_{\partial\Omega_{0}^{N}}\mathbf{T}\cdot\varphi
d\mathbf{S}.$ (3)

Above, $\Omega_{0}\subset R^{d}$ (d=2,3) is the initial domain in
consideration, $W=W(\mathbf{F}_{e})$ is a strain-energy density function,
$J_{g}=\det{\mathbf{G}}>0$ is the local volumetric variation due to growth –
the “Jacobian” from the Lagrangian coordinate $\mathbf{X}$ to a stress-free
grown manifold $\varphi_{g}(\mathbf{X},t)$, $\mathbf{T}$ is the standard
boundary traction, and $\mathbf{N}$ is the outward normal on free boundary
$\partial\Omega_{0}^{N}$. $\nabla_{0}$ is the gradient in the Lagrangian
coordinate, $\bm{\Pi}$ is called the 1st Piola–Kirchhoff stress tensor,
related to the Cauchy stress by $\bm{\sigma}=J^{-1}\bm{\Pi}\mathbf{F}^{T}$
with $J=\det\mathbf{F}>0$. In summary, the evolution of
$\mathbf{G}(\mathbf{X},t)$ drives the dynamics of $\varphi(\mathbf{X},t)$
through the minimization of Eq.(3) at each time.

As there have been a scarcity in the study on solution existence of the theory
in a general setting (see section 2.6 in [1]), most previous works derived
numerical and perturbative analytical methods to treat Eq.(2) case by case– in
a combination of different $\mathbf{G}(\mathbf{X},t)$, energy density function
$W(\mathbf{F}_{e})$, geometry of $\Omega_{0}$, and the associated boundary
conditions. In the case of perturbative analysis, they often treat a domain
with simple geometry (e.g., circular disk, annulus, spherical shell, etc) and
incompressibility constraint, which assist to find a trivial symmetric
solution of $\varphi$ to perform a buckling linear stability analysis (LSA)
(e.g.,[9, 5, 17, 15, 16, 7, 22]). A bifurcation threshold of
$\mathbf{G}(\mathbf{X},t)$ is found when the LSA shows that the trivial
solution becomes unstable. When the bifurcation is supercritical, the LSA can
also provide the selected wavenumber of the buckling.

The studies of the nonlinear solutions in the post-buckling stages (e.g.,[15,
16, 12, 18, 3]), and those starting with non-trivial geometry (e.g., smooth
brain before cortical folding [6, 21]), are mainly approached by finite
element methods (FEMs). While they have provide sophisticated simulations
which matches the experiments excitingly well with specific parameter
settings, they have not yet provided theoretical insights on these
simulations. For example, when smooth wrinkles develop to sharp creases, it is
still not clear if the two patterns correspond to the same solution, or there
has been transitions between multiple solutions.

To improve this, an eigenvalue analysis of the stiffness matrix of FEM has
been introduced by [11] to identify the bifurcation point of the shape
transitions. The eigen-directions of negative eigenvalues can identify the
instability modes, so that FEM does not have to be paired with LSA to initiate
the initial buckling. However, this approach can only improve our
understanding in supercritical transitions when the new solution in proximity
of the old solution can be approached guided by the eigen-directions followed
by Newton’s iterations.

In this study, we explore further the capacity of FEM and FEM-based
bifurcation analysis on growth-elasticity models in 2D. Besides using the
commercial software Abaqus FEA (e.g.,[15, 16, 12, 3]), previous FEMs solve the
mechanical equilibrium Eq.(2) by Galerkin approximation followed by Newton’s
iterations (e.g.,[18]) or by the discretizing the energy Eq.(3) followed by a
residual-force damping algorithm [20, 21]. Using the residual-force damping
algorithm alone lose insights on the solution stability and transitions due to
the lack of information on the stiffness matrix, while Newton’s iteration is
only reliable for local solution search. Here, we design a method combining
the strength of previous methods with further exploitation of the Hessian
matrix. By minimizing the energy Eq.(3) after uniform-strain discretization,
we use Hessian matrix of the energy (equivalent to stiffness matrix) not only
to guide the perturbation of a saddle solution to find new local minimizers
nearby, but also design an energy-stable gradient-descent algorithm to force
the solution to go to a desired groove when the local energy landscape is
rocky with multiple saddle points and minimizers. The solution is further
refined by Newton’s iterations. We demonstrate the combination of this energy-
stable algorithm with Newton’s iteration enables us to find multiple solutions
robustly in different types of bifurcations, especially with subcritical
bifurcations, or when multiple solutions are very close to each other.

The paper is organized as follows. In Sec. 2, we present a growth-elasticity
model with neo-Hookean elasticity which we aim to solve. In Sec. 3, we present
the numerical method with its discrete formulation and the solution search
algorithms. In Sec. 4, we apply the numerical method to study the shape
transitions in a bilayer annulus where the interior layer grows while being
constrained by the exterior layer. By adjusting the modulus ratio between the
two layers, we identify a variety of non-trivial solutions, including the
periodical smooth wrinkle-shape solutions born from supercritical bifurcations
[18], non-smooth crease solutions born subcritically [12], a transitory
aperiodical smooth solution due to the exchange of local minimality between
two close periodical wrinkle solutions, and aperiodical crease patterns before
periodical patterns occur. The latter two solutions have not been reported
previously. We conclude in Sec. 5 with discussions with previous works in [15,
18, 3, 12].

##  2 The model

In this study, we choose a neo-Hookean strain energy used in [21]

$\displaystyle
W(\mathbf{F}_{e})={\mu}/{2}(I_{3}^{-2/d}I_{1}-d)+{K}/{2}(I_{3}-1)^{2}$ (4)

for its simplicity. It involves two strain invariants
$I_{1}=\operatorname{Tr}(\mathbf{F}_{e}^{T}\mathbf{F}_{e})$ and
$I_{3}=\det\mathbf{F}_{e}$ with shear modulus $\mu$ and bulk modulus $K$,
where $d$ is the dimension. We note the strict incompressible neo-Hookean
energy has been frequently considered [9, 5, 17, 15, 16, 7, 22]:
$W(\mathbf{F}_{e})={\mu}/{2}(I_{1}-d)-p(I_{3}-1)$ with the Lagrangian-
multiplier pressure $p$ to ensure incompressibility $I_{3}=1$. We focus on the
compressible case in this study due to the extra challenge to impose the
constraint $I_{3}=1$ during optimization. Without the strict-incompressibility
constraint, the numerics in this study can be easily extended to consider more
complicated strain energy in different forms of $I_{1}$ and $I_{3}$ or
including other strain invariants [10].

We summarize the model as a variational problem. For any $0\leq t\leq
T_{end}$, given a trajectory of $\mathbf{G}(\mathbf{X},t)$ for
$\mathbf{X}\in\Omega_{0}$ and $\mathbf{T}(\mathbf{X},t)$ on the free boundary
for $\mathbf{X}\in\partial\Omega_{0}^{D}$, the field $\varphi(\mathbf{X},t)$
is a minimizer of Eq.(3) with

$\displaystyle\mathbf{F}_{e}=\nabla_{0}\varphi\mathbf{G}^{-1},$ (5)

the energy density function $W(\mathbf{F}_{e})$ from Eq.(4), and

$\displaystyle\varphi(\mathbf{X},t)=\mathbf{x}_{0}^{D}$ (6)

on the fixed boundary for $\mathbf{X}\in\partial\Omega_{0}^{D}$.

Given the above formulation, the dependence on time is merely through growth
$\mathbf{G}(\mathbf{X},t)$ and traction $\mathbf{T}(\mathbf{X},t)$, which
means at each time point $t_{j}\in[0,T_{end}]$, we solve
$\varphi(\mathbf{X},t_{j})$ from minimizing Eq.(3), using
$\mathbf{G}(\mathbf{X},t_{j})$ and $\mathbf{T}(\mathbf{X},t_{j})$.

##  3 Method Formulation

Figure 1: $\varphi^{k}\to\varphi^{k-,\alpha_{k}}\to\varphi^{k+,\alpha_{k}}$
are in the counterclockwise order in each triangle $\alpha$.

At each $t_{j}$, we consider the approximation of $\varphi$ in the space of
linear continuous interpolation functions, per triangulation of the initial
domain $\Omega(t_{0})=\cup_{\alpha}\Omega^{\alpha}(t_{0})\in R^{2}$. Thus, 1)
the distribution of $\varphi(t_{j})$ is determined by
$\\{\varphi^{k}(t_{j})\\}$ where $k$ indexes all the nodes
$\varphi^{k}(t_{j})=(\varphi^{k}_{x}(t_{j}),\varphi^{k}_{y}(t_{j}))\in\Omega(t_{j})$;
and 2) the deformation gradient
$\mathbf{F}^{\alpha}=[\nabla_{0}{\varphi}]^{\alpha}$ is uniform in each
$\Omega^{\alpha}$, determined by the coordinates of its vertices. In Fig. 1,
$\mathbf{F}^{\alpha}=\mathbf{m}^{\alpha}[\mathbf{M}^{\alpha}]^{-1}$ where the
two matrices
$\mathbf{m}^{\alpha}=[\varphi^{k-,\alpha}(t_{j})-\varphi^{k}(t_{j}),\varphi^{k+,\alpha}(t_{j})-\varphi^{k}(t_{j})]$
and
$\mathbf{M}^{\alpha}=[\varphi^{k-,\alpha}(t_{0})-\varphi^{k}(t_{0}),\varphi^{k+,\alpha}(t_{0})-\varphi^{k}(t_{0})]$.

For each $\mathbf{F}^{\alpha}$, we maintain the list of the three nodes in the
counterclockwise direction, such that $\det\mathbf{M}^{\alpha}>0$ initially.
Notice $\det\mathbf{M}^{\alpha}=2A^{\alpha}(t_{0}):=2\times\text{area of
}\Omega^{\alpha}(t_{0})$. At each time $t_{j}>t_{0}$, the sign of
$\det\mathbf{m}^{\alpha}=2A^{\alpha}(t_{j})$ informs us if
$\det\mathbf{F}^{\alpha}>0$ stay true. Within each triangle $\Omega^{\alpha}$,
we define the average elastic deformation
$\mathbf{F}_{e}^{\alpha}=\mathbf{F}^{\alpha}[\mathbf{G}^{\alpha}]^{-1}$ given
$\mathbf{G}^{\alpha}$ being the averaged $\mathbf{G}(\mathbf{X},t_{j})$ in
$\Omega^{\alpha}$.

###  3.1 Uniform-strain Energy Discretization

Given the triangulation and uniform-strain setting, we discretize the energy
functional by

$\displaystyle E$ $\displaystyle=$
$\displaystyle\sum_{\alpha}E^{\alpha}+\sum_{\beta}E^{\beta}$

where $\sum_{\alpha}E^{\alpha}$ accounts for
$\int_{\Omega_{0}}J_{g}W(\mathbf{F}_{e})d\mathbf{X}$ in the bulk and
$\sum_{\beta}E^{\beta}$ accounts for
$-\int_{\partial\Omega^{N}_{0}}\mathbf{T}\cdot\varphi d\mathbf{S}$ along the
free boundary ($\partial\Omega_{0}^{N}=\cup_{\beta}\partial\Omega^{\beta}$
where each $\partial\Omega^{\beta}$ are linear segments). For each triangle
$\Omega^{\alpha}$,

$\displaystyle
E^{\alpha}(\varphi^{1,\alpha},\varphi^{2,\alpha},\varphi^{3,\alpha},\mathbf{G}^{\alpha}):=J_{g}^{\alpha}W^{\alpha}A_{0}^{\alpha}$

where $J_{g}^{\alpha}=\det\mathbf{G}^{\alpha}$,
$W^{\alpha}=W(\mathbf{F}_{e}^{\alpha})=W(\mathbf{m}^{\alpha}[\mathbf{M}^{\alpha}]^{-1}[\mathbf{G}^{\alpha}]^{-1})$,
and $A^{\alpha}=\frac{1}{2}\det\mathbf{M}^{\alpha}$; and for the free-boundary
linear segments $\partial\Omega^{\beta}$ with two nodes
$\varphi^{1,\beta},\varphi^{2,\beta}$ and reference length $l_{0}^{\beta}$,

$\displaystyle
E^{\beta}(\varphi^{1,\beta},\varphi^{2,\beta}):=\frac{1}{2}\mathbf{T}^{\beta}\cdot(\varphi^{1,\beta}+\varphi^{2,\beta})l_{0}^{\beta}$

where $\mathbf{T}^{\alpha}$ is the averaged $\mathbf{T}({\mathbf{X},t_{j}})$
on $\partial\Omega^{\beta}$. For nodes on the fixed boundary
$\varphi^{k}(t_{0})\in\partial\Omega_{0}^{D}$, we have
$\varphi^{k}(t_{j})=\varphi^{k}(t_{0})$. To summarize, we write
$E(\varphi^{1},\varphi^{2},...,\varphi^{N})$ as a multivariable function of
the coordinates of the free nodes. From now on, we define
$\bm{\varphi}=[\varphi^{1},\varphi^{2},...,\varphi^{N}]$ concatenating the
$x-$ and $y-$ components of free $\varphi^{k}$’s.

###  3.2 Gradient Components

A local minimizer of $E(\bm{\varphi})$ must satisfy $\partial
E/\partial\bm{\varphi}=\mathbf{0}$. Here we provide the analytical formula
$\frac{\partial E}{\partial\varphi^{k}}\in R^{2}$ on each node. By denoting
$\alpha_{k}$ and $\beta_{k}$ for the triangles and boundary linear segment
connected to node $k$, respectively, we have

$\displaystyle\frac{\partial E}{\partial\varphi^{k}}$ $\displaystyle=$
$\displaystyle\sum_{\alpha_{k}}\frac{\partial
E^{\alpha_{k}}}{\partial\varphi^{k}}+\sum_{\beta_{k}}\frac{\partial
E^{\beta_{k}}}{\partial\varphi^{k}}$ (7)

where

$\displaystyle\frac{\partial E^{\beta_{k}}}{\partial\varphi^{k}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\mathbf{T}^{\beta_{k}}l_{0}^{\beta_{k}},$ (8)
$\displaystyle\frac{\partial E^{\alpha_{k}}}{\partial\varphi^{k}}$
$\displaystyle=$ $\displaystyle(J_{g}^{\alpha_{k}}\underbrace{\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}[\mathbf{F}_{e}^{\alpha_{k}}]^{T}}_{J_{e}^{\alpha_{k}}\bm{\sigma}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1})A_{0}^{\alpha_{k}}$
(9) $\displaystyle=$
$\displaystyle\frac{1}{2}\bm{\sigma}^{\alpha_{k}}\varphi^{k\perp,\alpha_{k}}.$

$\varphi^{k\perp,\alpha_{k}}$ is the counterclockwise orthogonal rotation of
$(\varphi^{k+,\alpha_{k}}-\varphi^{k-,\alpha_{k}})$ in the triangle
$\alpha_{k}$. See more details of the derivation steps in the Appendix. We
note this equation is equivalent to the two residual equations on each node
from the Bubnov–Galerkin method using linear continuous basis/interpolation
functions. Interestingly, this formula is more straightforward to interpret
the local discretized mechanics, since the force on each node $-\frac{\partial
E^{\alpha_{k}}}{\partial\varphi_{k}}=-\frac{1}{2}\bm{\sigma}^{\alpha_{k}}\varphi^{k\perp,\alpha_{k}}=\bm{\sigma}^{\alpha_{k}}\mathbf{n}(\frac{1}{2}|\varphi^{k\perp,\alpha_{k}}|)$
where $\mathbf{n}=-\varphi^{k\perp,\alpha_{k}}/|\varphi^{k\perp,\alpha_{k}}|$
is the outward unit normal on the edge
$(\varphi^{k+,\alpha_{k}}-\varphi^{k-,\alpha_{k}})$.

###  3.3 Hessian Components

Besides the condition $\partial E/\partial\bm{\varphi}=\mathbf{0}$, a local
minimizer of $E(\bm{\varphi})$ must have its Hessian matrix
$\partial^{2}E/(\partial\bm{\varphi}\partial\bm{\varphi})$ be positive semi-
definite. The Hessian matrix is equivalent to the stiffness matrix in previous
Galerkin-FEMs, which seemed to be computed by symbolic operations or numerical
differentiation of the residual equations, since they have not been provided
explicitly to the best of our knowledge. Here we provide the analytics of
Hessian components between two nodes $k_{1}$ and $k_{2}$ via triangle
$\alpha$, which we use to assemble
$\partial^{2}E/(\partial\bm{\varphi}\partial\bm{\varphi})$:

$\displaystyle\frac{\partial^{2}E^{\alpha}}{\partial\varphi^{k_{1}}\partial\varphi^{k_{2}}}=\frac{\mu}{2}\bigg{(}\frac{\partial^{2}E_{\mu}^{\alpha}}{\partial\varphi^{k_{1}}\partial\varphi^{k_{2}}}\bigg{)}+\frac{K}{2}\bigg{(}\frac{\partial^{2}E^{\alpha}_{K}}{\partial\varphi^{k_{1}}\partial\varphi^{k_{2}}}\bigg{)}$
(10)

where

$\displaystyle\frac{\partial^{2}E^{\alpha}_{\mu}}{\partial\varphi^{k_{1}}\partial\varphi^{k_{2}}}=\frac{(\varphi^{k_{1}\perp}\cdot\varphi^{k_{2}\perp})}{\det[\mathbf{m}^{\alpha}]}[\mathbf{B}^{\alpha}]^{-1}-\frac{\text{tr}[\mathbf{B}^{\alpha}]\epsilon(\varphi^{k_{1}},\varphi^{k_{2}})}{2\det{[\mathbf{B}^{\alpha}]}}\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}$ (11)

and

$\displaystyle\frac{\partial^{2}E^{\alpha}_{K}}{\partial\varphi^{k_{1}}\partial\varphi^{k_{2}}}=\frac{1}{J_{g}^{\alpha}\det[\mathbf{M}^{\alpha}]}[\varphi^{k_{1}\perp}\otimes\varphi^{k_{2}\perp}]+(J_{e}^{\alpha}-1)\epsilon(\varphi^{k_{1}},\varphi^{k_{2}})\begin{bmatrix}0&-1\\\
1&0\end{bmatrix}$ (12)

where
$[\mathbf{B}^{\alpha}]=[\mathbf{F}_{e}^{\alpha}][\mathbf{F}_{e}^{\alpha}]^{T}$
and $\epsilon(\varphi^{k_{1}},\varphi^{k_{2}})=\left\\{\begin{array}[]{ c l
}0&\quad\textrm{if }\varphi^{k_{1}}=\varphi^{k_{2}}\\\ -1&\quad\textrm{if
}\varphi^{k_{1}}=\varphi^{k_{2}+,\alpha}\\\ 1&\quad\textrm{if
}\varphi^{k_{1}}=\varphi^{k_{2}-,\alpha}\end{array}.\right.$

###  3.4 Iterative Optimizations

With the components of both the gradient $\partial E/\partial\bm{\varphi}$ and
the Hessian $\partial^{2}E/(\partial\bm{\varphi}\partial\bm{\varphi})$ being
formulated, we are ready to utilize optimization algorithms to search for
local minimizers. In this study, we have designed an algorithm to utilize
Newton’s method and a customized energy-stable gradient-descent algorithm,
where the former is used to accelerate convergence when the numerical iterates
are close to a critical point, and the latter is used to provide initial guess
to the Newton’s iterates. The energy-stable gradient-descent steps are
necessary for the robustness of whole algorithm close to a bifurcation, as we
will demonstrate in the Results. Newton’s method gives the iteration:

$\displaystyle\bm{\varphi}_{i+1}-\bm{\varphi}_{i}=-[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]^{-1}\partial
E/\partial\bm{\varphi}_{i}.$ (13)

The gradient-descent algorithm gives the iteration:

$\displaystyle\bm{\varphi}_{i+1}-\bm{\varphi}_{i}=-\Delta s\partial
E/\partial\bm{\varphi}_{i}.$ (14)

As long as the maximal Hessian-eigenvalue
${\lambda_{max}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]}>0$
($\bm{\varphi}_{i}$ is not a local maximizer), we choose $\Delta s$ such that
the energy is non-increasing:

$\displaystyle\Delta s\leq\Delta
s_{stable}:={2}/{\lambda_{max}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]}.$
(15)

See the derivation of $\Delta s_{stable}$ in the Appendix. Notice if
$\bm{\varphi}_{i}$ is a local maximizer, any small $\Delta s>0$ will suffice
(also see Appendix).

Acknowledging the computational expense of this search, it nonetheless helps
us to localize the Newton’s iteration in a region of interest. This strategy
is especially important when we are searching for a specific solution among
two or multiple solutions. To improve the efficiency, we also follow the step
size

$\displaystyle\Delta{s}_{opt}={2}/\bigg{(}\lambda_{max}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]+\lambda_{min}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]\bigg{)}$
(16)

to achieve the optimal 1st-order rate of convergence

$\lambda_{opt}=\frac{\lambda_{max}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]-\lambda_{min}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]}{\lambda_{max}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]+\lambda_{min}[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]}\leq
1$ (17)

when $[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]\geq
0$ is positive semi-definite. The derivation of $\Delta{s}_{opt}$ is also
provided in the appendix. At each gradient-descent iteration, we choose
$\Delta{s}=\min{(\Delta{s}_{opt},\Delta s_{stable})}$. The solution search
protocol is summarized in Algorithm 1.

Data: $[\mathbf{G}]^{\alpha}(t_{j})$ and $\mathbf{T}^{\beta}(t_{j})$.

Result: $\bm{\varphi}(t_{j})$.

Start with the initial guess $\bm{\varphi}_{0}$;

while  _$||\partial E/\partial\bm{\varphi}_{i}||_{2} >Tol$ and
$Iter<MaxIter_{1}$_ do

$\bm{\varphi}_{i+1}=\bm{\varphi}_{i}-\Delta s\partial
E/\partial\bm{\varphi}_{i}$ with $\Delta{s}=\min{(\Delta{s}_{opt},\Delta
s_{stable})}$;

end while

while  _$||\partial E/\partial\bm{\varphi}_{i}||_{2} >Tol$ and
$Iter<MaxIter_{2}$_ do

$\bm{\varphi}_{i+1}=\bm{\varphi}_{i}-[\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})]^{-1}\partial
E/\partial\bm{\varphi}_{i}.$

end while

Algorithm 1 An iterative algorithm at each time $t_{j}$. See texts for the
choice of initial guess $\bm{\varphi}_{0}$, $Tol$, $MaxIter_{1}$, and
$MaxIter_{2}$.

In this study, we choose the error tolerance to be $Tol=1e-7$. In most
results, the numerical solutions meet the error tolerance under
$MaxIter_{1}=1e4$ and $MaxIter_{2}=20$. For particularly challenging cases, we
need to adjust $MaxIter_{1}=1e5-1e6$ to locate an initial guess for the
Newton’s iterates. The selection of $\bm{\varphi}_{0}$ at each time is either
based on the solution from the last time step or a perturbed guess utilizing
eigenpair information from the Hessian. Further details will be provided in
Results.

##  4 Results

Figure 2: Schematics of the initial domain.

We apply Algorithm 1 to study shape transitions in a bilayer annulus where the
interior layer grows isotropically and uniformly while being constrained by
the exterior layer without any growth (see Fig. 2). The interior boundary is
free from tractions ($\mathbf{T}=\mathbf{0}$) and the exterior boundary is
fixed ($\varphi(\mathbf{X},t)=\mathbf{X}$). This set up is relatively simple,
but the corresponding numerical solutions present a rich set of behaviors
–wrinkles, smooth folds, creases [15, 12, 18, 3] (more discussions in the
Conclusion).

We consider an annular region with its initial inner radius half of its outer
radius, $R_{in}=0.5$ and $R_{out}=1$, and triangulate it by $N_{x}=12$ layers
of triangles along the radius, and $N_{y}=92$ layers of triangles along the
circumference. These triangles are similar, and increase their size from the
interior boundary to the exterior boundary by a ratio of
$(R_{out}/R_{in})^{1/N_{x}}\approx 1.06$ between successive layers. We set the
growing domain to be the first two interior layer of triangles (see Fig.3a).
Since the growth is isotropic, $\mathbf{G}=g\mathbf{I}_{2}$, we have a single
parameter $g$ to control the elastostatic system. As such, we study the system
as $g$ increases from $g=1$ (no growth) in the growing domain, in three
different cases with different elastic-modulus ratio between the growing and
non-growing domain (i.e., varying the ratio between the shear modulus
$\mu_{g}$ in the growing domain and $\mu_{ng}$ in the non-growing domain; and
the ratio between the bulk modulus $K_{g}$ in the growing domain and $K_{ng}$
in the non-growing domain). We note that this study is not about
systematically studying the effect of $R_{out}/R_{in}$ in combination of
$\mu_{ng}/\mu_{g}$ and $K_{ng}/K_{g}$, but to demonstrate the capacity of the
new method to produce new solution details in various scenarios.

Figure 3: (a) Three snapshots of the evolution of the annular solution
$s^{0}$. “s” at the center of a configuration indicates it is a saddle point,
not a local minimizer. (b) The three smallest eigenvalues of the Hessian of
$s^{0}$ versus $g$. The vector fields are generated by plotting the
eigenvector components on the corresponding nodes.

###  4.1 Case 1: $\mu_{ng}/\mu_{g}=0.1$ and $K_{ng}/K_{g}=0.1$

As $g$ increases from unity, the annular solution $\bm{\varphi}(s^{0})$ soon
becomes a saddle point at around $g=1.261$ (Fig.3b). Close to the transition,
the smallest eigenvalue of the Hessian matrix $H(s^{0})$ is associated with an
eigenplane $E_{\lambda_{min}}[H(s^{0})]=N(H(s^{0})-\lambda_{min}\mathbf{I})$,
where $N(\Box)$ denotes the nullspace of $\Box$. Any vector $v_{5}^{0}\in
E_{\lambda_{min}}[H(s^{0})]$ delineates a vector field of a wavenumber $m=5$
around the circumference of the annular domain (see the blue vector field in
Fig.3b).

The two-dimensional $E_{\lambda_{min}}[H(s^{0})]$ indicates a family of
shallow-wrinkled solutions with a wavenumber $m=5$. We compute one of the
solutions by perturbing $\bm{\varphi}(s^{0})$ in the direction of a
$v_{5}^{0}\in E_{\lambda_{min}}[H(s^{0})]$ with a small magnitude, followed by
Algorithm 1. We show five samples of these solutions $\\{s^{5}\\}$ at
$g=1.262$ with a rotation of $72^{\circ}$ clockwisely (Fig.4a, first row)
followed by their development at $g=1.284$ (Fig.4a, second row). At each value
of $g$, the family of solutions $\\{s^{5}\\}$ have equal distance
$d(g)=||\varphi(s^{5})-\varphi(s^{0})||_{2}$ from $s^{0}$, which increases as
$g$ increases. By fitting the $g$ versus $d$ by quadratic polynomials, we find
the first-order coefficient is negligible compared to the second-order
coefficient (see Fig.4b), suggesting a relation $(g-1.261)\propto d^{2}$. This
indicates that the family $\\{s^{5}\\}$ is born from a pitchfork-like
supercritical bifurcation at around $g_{c}=1.261$ , and this relation $d$
versus $g$ is analogous to the relation $\epsilon\propto\sqrt{g-g_{c}}$ in the
weakly nonlinear analysis where $\epsilon$ is the buckling amplitude [2, 4].

Figure 4: (a) Five isoenergetic solutions from $\\{s^{5}\\}$ at $g=1.262$ and
$g=1.284$. (b) The growth stretch parameter $g$ versus
$d=||\bm{\varphi}(s^{5})-\bm{\varphi}(s^{0})||_{2}$. (c) The force-residual
$||\partial E/\partial\bm{\varphi}_{linearized}||_{2}$ due to the linear
approximation versus $g$; and the distance between the nonlinear and linear
solutions $\sqrt{d^{2}-d_{H}^{2}}/d$ relative to the distance between the
nonlinear solution $\\{s^{5}\\}$ and $s^{0}$.

To measure the nonlinearity of these solutions (how inaccurately they can be
approximated by the linear approximation), we first identify the linear
approximations as the projection of the full nonlinear solutions on the
eigenplane $E_{\lambda_{min}}[H(s^{0})]$. Due to the equal distance from
$s^{0}$, the projection is a circle with radius $d_{H}$ in the eigenplane
$E_{\lambda_{min}}[H(s^{0})]$ centered at $s^{0}$. We quantify the
approximation error by 1) the relative distance $\sqrt{d^{2}-d_{H}^{2}}/d$ and
2) the gradient error $||\partial E/\partial\bm{\varphi}_{linearized}||_{2}$
from the linear approximation. As expected, both increases as $g$ increases
(Fig. 4c).

Figure 5: (a,b) The configurations of a $\\{s^{6}\\}$ and $\\{s^{4}\\}$
solution at $g=1.286$, respectively, juxtaposed with the vector fields
corresponding to the smallest Hessian-eigenvalues. “s” at the center of a
configuration indicates it is a saddle point, not a local minimizer. Both
$\\{s^{6}\\}$ and $\\{s^{4}\\}$ have their smallest eigenvalue negative. Each
both case, two different orthogonal vector fields (blue and orange) from the
eigenplane are plotted. They all resemble an undulation of $m=5$. (c) The
configuration of a $\\{s^{5}\\}$ juxtaposed with the smallest-eigenvalue
vector field, which is nearly zero. This vector field depicts the local
tangent of $\\{s^{5}\\}$-solutions at the current solution $s^{5}$.

As $g$ further increases, more eigenvalues of $H(s^{0})$ become negative, and
they are also associated with two-dimensional eigenplanes. The vector fields
generated by eigenvectors $v_{4}^{0}$ and $v_{6}^{0}$ have wavenumbers of
$m=6$ and $m=4$, respectively, along the circumference of $s^{0}$ (see the
green and red vector fields in Fig.3b). We find the corresponding solutions
$\\{s^{6}\\}$ and $\\{s^{4}\\}$ by performing the eigen-direction-guided
perturbations followed by Algorithm 1, except these solutions turn out to be
saddle points at birth. This makes sense if $\\{s^{6}\\}$ and $\\{s^{4}\\}$
are born from $s^{0}$ when $s^{5}$ is still close enough to $s^{0}$. We
confirm this proximity by checking the negative-value eigenspace of $H(s^{6})$
and $H(s^{4})$, and we find both vector fields echo $\\{s^{5}\\}$ (Fig.5a and
b). In contrast, the smallest eigenvalue of the locally stable $\\{s^{5}\\}$
is vanishing (see Fig.5c). Its associated eigenvector approximates the kernel
of $H(s^{5})$, which is the local tangent of $\\{s^{5}\\}$.

Figure 6: (a) The evolution of a $\\{s^{5}\\}$-, $\\{s^{6}\\}$-, and
$\\{s^{4}\\}$-solution versus $g$. (b) The smallest eigenvalue of the
$\\{s^{5}\\}$-, $\\{s^{6}\\}$-, and $\\{s^{4}\\}$-solution versus $g$. (c) The
energy of the $\\{s^{5}\\}$-, $\\{s^{6}\\}$-, and $\\{s^{4}\\}$-solution
relative to $E(s^{0})$ versus $g$. “s” at the center of a configuration
indicates it is a saddle point, not a local minimizer.

We find that $\\{s^{6}\\}$ and $\\{s^{4}\\}$ both become locally stable as $g$
further increase (Fig.6a and b). However, $\\{s^{5}\\}$ is still the absolute
minimizer among the three family of solutions up to $g=2$ (Fig.6c). Strictly
speaking, we cannot exclude the existence of other local minimizers that have
a lower energy than $\\{s^{5}\\}$. We do not further increase $g$, as self-
contact occur both in $\\{s^{5}\\}$ and $\\{s^{4}\\}$. We leave the
consideration of treating self contact [21] in future work. We also remark
that the smallest eigenvalue of $\\{s^{5}\\}$ becomes positive when $g$
becomes $1.6$. This indicates the triangulation cannot accurately represent
smooth solutions anymore, and mesh refinement is needed. We also leave this to
future works for post-buckling computations.

###  4.2 Case 2: $\mu_{ng}/\mu_{g}=0.25$ and $K_{ng}/K_{g}=0.25$

In this case, we find the annular solution $\bm{\varphi}(s^{0})$ becomes a
saddle point at $g=1.475$ by increasing $g$ from 1 with $\Delta g=1e-3$. The
same as Case 1, the negative eigenvalue at $g=1.475$ is associated with a
2-dimensional eigenplane, whose vector fields present a wavenumber of $m=5$.
Following the same strategy as Case 1, perturbing $s^{0}$ by negative eigen-
direction and applying Algorithm 1, we find the first nontrivial solution
$\\{s^{5}\\}$ (Fig.7a). Immediately at $g=1.476$, $H(s^{0})$ has its second
smallest eigenvalue becoming negative, which results in the second nontrivial
solution $\\{s^{6}\\}$. The same as Case 1, $\\{s^{6}\\}$ is born as a saddle
point (Fig.7a,d). Both solutions are born from a subcritical bifurcation
(Fig.7b,c). In contrast to Case 1, the two new solutions are born very closely
to each other.

Figure 7: (a) The evolution of a $\\{s^{5}\\}$-, $\\{s^{6}\\}$-, and
$\\{s_{1}\\}$-solution versus $g$. Red boxes highlight the solution with the
smallest energy at each value of $g$. “s” at the center of a configuration
indicates it is a saddle point, not a local minimizer. The purple vector field
on $s^{6}(g=1.494)$ is from the eigenvector associated with the negative
eigenvalue of $H(s^{6})$. (b) The energy of the $\\{s^{5}\\}$-,
$\\{s^{6}\\}$-, and $\\{s_{1}\\}$-solution relative to $E(s^{0})$ versus $g$.
(c) The distance from the annular solution
$||\bm{\varphi}(s^{5})-\bm{\varphi}(s^{0})||_{2}$,
$||\bm{\varphi}(s^{6})-\bm{\varphi}(s^{0})||_{2}$, and
$||\bm{\varphi}(s_{1})-\bm{\varphi}(s^{0})||_{2}$ versus $g$. (d) The smallest
eigenvalue for $H(s^{5})$, $H(s^{6})$, and $H(s_{1})$, respectively, versus
$g$.

As $g$ further increases with $\Delta g=1e-3$, $\\{s^{5}\\}$ and $\\{s^{6}\\}$
switches stability at around $g=1.477$. At $g=1.477$, both solutions are
saddle points. By perturbing either $\\{s^{5}\\}$ or $\\{s^{6}\\}$ by vector
fields from their unstable eigenplanes followed by Algorithm 1, we obtain the
same new solution $\\{s_{*}\\}$ (Fig.8a). We emphasize that since
$\\{s^{5}\\}$ or $\\{s^{6}\\}$ are very close already, using Newton’s method
alone fails to identify $\\{s_{*}\\}$–when we started with one solution, say
$s^{5}$, we end up finding the other, in this case $s^{6}$. The energy-stable
gradient-descent steps are necessary to localize the search to the valley
between $\\{s^{5}\\}$ and $\\{s^{6}\\}$.

Figure 8: (a) The evolution of a $\\{s_{*}\\}$-solution versus $g$,
transitioning from a shape of $\\{s^{5}\\}$ to $\\{s^{6}\\}$. (b) The smallest
eigenvalue for $H(s_{*})$ versus $g$, juxtaposed with those of $H(s^{5})$ and
$H(s^{6})$. $\\{s_{*}\\}$ is locally stable while the other two are not. (c)
The distance of $\\{s_{*}\\}$ from the annular solution
$||\bm{\varphi}(s_{*})-\bm{\varphi}(s^{0})||_{2}$, juxtaposed with those of
$\\{s^{5}\\}$ to $\\{s^{6}\\}$. Notice
$||\bm{\varphi}(s_{*})-\bm{\varphi}(s^{0})||_{2}$ converges to
$||\bm{\varphi}(s^{5})-\bm{\varphi}(s^{0})||_{2}$ at the lower end of $g$, and
converges to $||\bm{\varphi}(s^{6})-\bm{\varphi}(s^{0})||_{2}$ at the higher
end. (d) The energy $E(s_{*})$ versus $g$, juxtaposed with those of $E(s^{5})$
and $E(s^{6})$. While $\\{s^{5}\\}$ and $\\{s^{6}\\}$ are exchanging the
minimality of energy among they two, $\\{s_{*}\\}$ has the lowest energy
throughout the exchange.

We further study $\\{s_{*}\\}$ around $g=1.477$ with $\Delta g=1e-4$. We find
that $\\{s_{*}\\}$ is a temporary solution that exists in an interval of $g$
within $[1.4765,1.4773]$. The shape transition in this window suggests that it
is a solution born from $s^{5}$ and later merges with $s^{6}$. We further
study its smallest eigenvalue $\lambda_{min}(s_{*})$ (Fig.8b), the distance
from $s^{0}$ (Fig.8c), and its energy relative to $s^{0}$ (Fig.8d), and
compare these values with counterparts in $\\{s^{5}\\}$ and $\\{s^{6}\\}$. All
these results agree with the suggestion that this transitory solution exits
only during the exchange of stability between $\\{s^{5}\\}$ or $\\{s^{6}\\}$.

As we continue to increase $g$ after when $\\{s^{6}\\}$ becomes stable, we
find $\\{s^{6}\\}$ becomes unstable again as $g=1.494$ (with an increment of
$\Delta g=1e-3$). Perturbing $\\{s^{6}\\}$ by the vector field from the
unstable eigenspace (see the purple vector field in Fig.7a) followed by
Algorithm 1, we find a new solution (family) $\\{s_{1}\\}$. Since the
configuration of $\\{s_{1}\\}$ resembles the vector field from the unstable
eigenspace, we have anticipated that this transition is supercritical from
$\\{s^{6}\\}$. To our surprise, as we decrease the value of $g$ on
$\\{s_{1}\\}$, we find it exists for smaller values of $g$ down to $g=1.478$.
This means $\\{s_{1}\\}$ coexists with $\\{s^{6}\\}$ as local minimizers
before $\\{s^{6}\\}$ becomes locally unstable. Thus the birth of $\\{s_{1}\\}$
is in fact subcritical from $\\{s^{6}\\}$. One may also question if
$\\{s_{1}\\}$ is born continuously from $\\{s^{0}\\}$. We can see from
Fig.7b,c that it is born with a finite distance from $s^{0}$ and has higher
energy than the trivial solution $s^{0}$ when it is born.

To catch the action that $\\{s_{1}\\}$ is born from deflecting the local
energy landscape at where it is born, we have done a refined search of $g$ in
$[1.47770725,1.47771725]$ with $\Delta g=5e-6$ and $[1.47772725,1.47779725]$
with $\Delta g=1e-5$. One can see the transition of $\lambda_{min}(s_{1})$
from $0$ is very sharp with respect to the change of $g$ (Fig.7d).

Similar to the smooth solutions such as $\\{s^{5}\\}$ and $\\{s^{6}\\}$, there
should be a family of solutions $\\{s_{1}\\}$ due to the energy invariance to
the rotational displacement of the creases, but we will not obtain a null
space for each $H(s_{1})$, since each solution in $\\{s_{1}\\}$ cannot feel
its neighbor – they are divided by a finite distance due to the size of the
finite elements. Different from the smooth wrinkle solutions, we anticipate
this property remains as the mesh refines, since there is a discontinuity in
$\nabla_{0}\varphi$ which cannot be resolved by mesh refinement. In general
for creases, they only have a non-empty null space in the Hessian at the
critical point when they are just born or when they become unstable again.

Figure 9: The evolution of crease formations from $\\{s_{1}\\}$ to
$\\{s_{6}\\}$. Red boxes highlight the solution with the smallest energy at
each value of $g$. “s” at the center of a configuration indicates it is a
saddle point, not a local minimizer. See inset for the energy versus $g$ for
each solution. Though not identical, the energy among two solutions at each
value of $g$ are close. In Case 3, the energy curves will be more splayed out.

As $g$ further increases, the smallest eigenvalue of $\\{s_{1}\\}$ decreases,
suggesting $\\{s_{1}\\}$ will become unstable at a larger critical value of
$g$ (Fig.7d). As we continue to increase $g$, we aim to find a solution path
of local minimizers rather than accurately identify all the bifurcation points
where one solution becomes unstable or when one solution is born. Following
this cost-efficient strategy, we find a $\\{s_{2}\\}$ from using the initial
configuration of $s_{1}(g=1.5)$ as $\bm{\varphi}_{0}$ in Algorithm 1 at
$g=1.51$ (Fig.9, row $s_{2}$). We emphasize again that the energy-stable
gradient-descent search is necessary to find $\\{s_{2}\\}$. We have tried to
use Newton’s method alone, and we fail to find $\\{s_{2}\\}$ since
$\\{s_{2}\\}$ is not local to $\\{s_{1}\\}$.

Then we decrease and increase $g$ with $\Delta g=5e-3$ on $s_{2}(g=1.51)$ and
find $\\{s_{2}\\}$ exists in the window $g\in[1.49,1.52]$. Similarly, we find
$s_{4}^{a}$ from using $s_{2}(g=1.52)$ as $\bm{\varphi}_{0}$ at $g=1.525$
(Fig.9, row $s_{4}^{a}$), and find $s_{6}$ from using $s_{4}^{a}(g=1.54)$ as
$\bm{\varphi}_{0}$ at $g=1.545$ (Fig.9, row $s_{6}$).

One can see from Fig.9a that each value of $g$ has at least two solutions,
among which we generate a minimal-energy path (solutions in red boxes). We
find more local minimizers such as $s^{5}(g=1.545)$ and $s^{7}(g=1.55)$, and a
saddle $s^{4}(g=1.55)$ as well as other aperiodical patterns with creases (not
shown). None of them has a lower energy than the minimal energy among
$\\{s_{1}\\}$, $\\{s_{2}\\}$, $\\{s_{4}^{a}\\}$, and $\\{s_{6}\\}$ at a
specific value of $g$ when they are available. Of course, we cannot exclude
the possibility that other minimizers exist with a lower energy. But at least
we have provided a path of shape transitions in the crease-forming phase of
the post-buckling evolution.

###  4.3 Case 3: $\mu_{ng}/\mu_{g}=1$ and $K_{ng}/K_{g}=1$

As we increase $g$ on the annular solution $s^{0}$ with $\Delta g=5e-3$, we
find it becomes a saddle at $g=1.7$ (see Fig.10, row $s_{0}$). The eigenplane
corresponding to the negative eigenvalue gives us a $m=6$ vector field. We
thus perturb $s^{0}$ according to the vector field, and find a local minimizer
$s_{6}$. Again, without the gradient-descent search, we cannot find this
solution.

Figure 10: The evolution of crease formations from $\\{s^{0}\\}$ to
$\\{s_{6}\\}$. Red boxes highlight the solution with the smallest energy at
each value of $g$. “s” at the center of a configuration indicates it is a
saddle point, not a local minimizer.

We will show in the following all the solutions we find are subcritical.
Similar to the crease-forming phase of the post-buckling study in Case 2, we
do not pursue the exact value of $g$ where the subcritical solutions occur or
become unstable. Instead, we focus on finding a path of local minimizers (with
the smallest energy among all the solutions we have found) to describe the
shape dynamics driven by the change of $g$.

Figure 11: (a) The relative energy $E/E(s_{0})$ versus $g$ among $s_{0}$,
$\\{s_{1}\\}$, $\\{s_{2}^{b}\\}$, $\\{s_{3}^{a}\\}$, $\\{s_{4}^{a}\\}$, and
$\\{s_{6}\\}$. (b) The distance of $s_{0}$, $\\{s_{1}\\}$, $\\{s_{2}^{b}\\}$,
$\\{s_{3}^{a}\\}$, $\\{s_{4}^{a}\\}$, and $\\{s_{6}\\}$ from $s_{0}$.

First, we decrease $g$ with $\Delta g=5e-3$ on $\\{s_{6}\\}$, and find
$\\{s_{6}\\}$ occurs as early as $g=1.65$ (see Fig.10, row $s_{6}$). At
$g=1.65$, close to when $\\{s_{6}\\}$ is born, we find it is away from $s_{0}$
by $d=||\varphi(s_{6})-\varphi(s_{0})||_{2}\approx 0.75$ with the energy
larger than $s_{0}$ (see Fig.11, a and b). These results confirm that it is
born from a subcritical bifurcation. Further decreasing $g$ on $s_{6}$, we
find a $s_{4}^{a}(g=1.645)$. By increasing and decreasing $g$ on $s_{4}^{a}$,
we find it is born near $g=1.63$ and continue to exist beyond $g=1.7$ (Fig.10,
row $s_{4}^{a}$). Similar to $s_{6}$, we show in Fig.11, a and b, that $s_{4}$
is born with $d\approx 0.6$ and $E(s_{4}^{a})/E(s_{0})>1$.

Suspecting that there exist other local minimizers with different number of
creases and their spatial arrangements within $1.61\leq g\leq 1.7$, we perturb
$s_{0}(g=1.7)$ with vector fields corresponding to the 2-7th smallest (yet
positive) eigenvalues of $H(s_{0})$ at $g=1.7$ with a finite magnitude
$|v|=1$. This magnitude is large enough to perturb the solution across an
energy barrier and locate a local minimizer in most cases – too small of a
perturbation only leads us back to $s_{0}$ using Algorithm 1.

The 2-5th smallest eigenvalue of $H(s_{0})$ is associated with a vector field
of a circumferential wave with wavenumber $m=\\{5,7,4,8\\}$, respectively. The
6th smallest eigenvalue corresponds to an angular-rotation vector field, and
the 7th smallest eigenvalue corresponds to a wave with $m=9$.

For a $m=4$ perturbation, we find a $s_{4}(g=1.7)$, which occurs as early as
$g=1.635$ (Fig.10, row $s_{4}$). For a $m=8$ perturbation, we still find a
$s_{4}(g=1.7)$. For the rotational perturbation, the solution goes back to
$s_{0}$.

We find more solutions from perturbations of $m=\\{5,7,9\\}$. A $m=5$
perturbation leads to a $s_{5}$, which occurs as early as $g=1.64$ (row
$s_{5}$). Further decreasing $g$ on $s_{5}(g=1.64)$ finds an aperiodical
$s_{3}^{a}$ (row $s_{3}^{a}$) that occurs as early as $g=1.625$. By decreasing
$g$ on $s_{3}^{a}(g=1.625)$, we find a $s_{1}$ that occurs as early as
$g=1.615$. In fact, $s_{1}$ is earliest solution we can find, and is closest
to $s_{0}$ (see Fig.11 a,b). Still, when it is born, it has a higher energy
than $s_{0}$.

A $m=7$ perturbation on $s_{0}(g=1.7)$ leads to a $s_{7}$, which occurs as
early as $g=1.675$. Further decreasing $g$ on $s_{7}(g=1.675)$ lands on an
aperiodical solution of $s_{4}^{b}$ which occurs as early as $g=1.64$. Further
decreasing $g$ on $s_{4}^{b}(g=1.64)$ finds us an aperiodical $s_{2}^{a}$. A
$m=9$ perturbation on $s_{0}(g=1.7)$ leads to an aperiodical $s_{4}^{c}$ which
occurs at around $g=1.675$. Further decreasing $g$ on $s_{4}^{c}(g=1.675)$
finds us an aperiodical $s_{3}^{b}$. Interestingly, further decreasing
$s_{3}^{b}(g=1.635)$ finds us a nearly periodical $s_{2}^{b}$ which occurs as
early as $g=1.62$ (row $s_{2}^{b}$ in Fig.10).

Choosing the absolute minimizer among these solutions at each value of $g$, we
find a path of the shape evolution where the number of crease increases by $1$
at each discontinuous shape transition, except from $s_{4}^{a}\to s_{6}$.
Again, we cannot claim to have exhausted all the solutions.

##  5 Conclusion

In this work, we have introduced an energy-optimization method to study
growth-elasticity models. This method obviously has its connection to previous
FEM methods, as mentioned in the introduction. The novelty and the strength of
this method come from its robustness in finding new solutions, owing to the
energy-stable gradient-descent algorithm preceding the Newton’s iterates. In a
long-range search in subcritical bifurcations, this algorithm does not blow up
when Hessian becomes near singular. In very localized searches when multiple
solutions are very close to each other, this algorithm forces the numerical
iterates to go down to the valley containing the desired solution and forbids
overshooting.

When the non-growing domain is softer than the growing domain (Case 1 and Case
2), we find the initial buckling as a family of smooth wrinkled solution,
which is similar to previous results from linear stability analysis ([17, 15],
though we have not considered a pressure at the interior surface) and FEMs
[18]. [18] has shown that the interior boundary develops smooth wrinkles which
later develops into finite-amplitude smooth folds with different wavenumbers.
However, the detail of the initial bifurcation is not captured by [18]. More
to [18], we fully reveal the supercritical nature of the initial transition,
and find that the smooth wrinkle solutions that occur after the initial
transition are born as saddle points, meaning they are not initially stable.
As growth continues, in Case 1, these saddle points can become locally stable,
but they do not become the absolute minimizer among the three solutions we
study.

In Case 2, when the first and second new solutions are very close to each
other, they two solutions switch stability as $g$ increases, which is also not
report before. Thanks to the energy-stable gradient-descent searches, we are
able to identify a new aperiodical minimizer between the two saddle solutions
during the switch of stability. The previous FEM [18] cannot catch such a
temporary minimizer since they are based on Newton’s method alone. We have
confirmed this by applying Newton’s method alone which always give us the
saddle solutions. To the best of our knowledge, such transitory aperiodical
solutions have not been reported in theory previously, especially when they
occur very close to the initial bifurcation. Notice this is different from
period-doubling patterns, which has been seen in the post-buckling stage in
previous work [15].

In Case 2, we find the wrinkled solutions become saddle points eventually thus
cannot further develop into smooth folds as in Case 1. Instead, we find
locally stable crease solutions. We confirm the subcritical nature of the
crease solutions by showing the next solution occurs before the previous
solution become unstable. For a $m=6$ smooth wrinkle solution to become a
periodical pattern with $6$ creases, we show multiple subcritical transitions
in between, which is different from previous simulations where all the
indentations due to smooth wrinkles develop into creases simultaneously [3].
However, there are multiple differences between [3] and this work, including
the 3D cylindrical geometry in [3] versus the 2D geometry in this work,
different strain-energy function, material property overlay, and the mesh
resolution. Future work will further clarify the nature of crease development
by mesh refinements.

In Case 3 when the non-growing and growing domains have the same elastic
moduli, we find the annular solution loses local stability after multiple non-
trivial solutions have been born from subcritical bifurcations. In this case,
a linear stability analysis approach would fail to predict the initial shape
transition. Previously the critical condition of crease formation has been
described by comparing the elastic energy in a prescribed creased solution and
that in a smooth solution [12]. When the two energies are equal, the crease
solution is considered to be selected. However, only a limited number of
periodical crease solutions have been found. Thanks to the robustness of our
method, by exploring the g-window for each periodical crease solutions, we
have harvested multiple aperiodical solutions close to the limit when these
periodical solutions become unavailable. Eventually we have paved a trajectory
of crease formation evolution as $g$ grows dynamically.

##  6 Appendix

###  6.1 The derivation of Eq.(9)

$\displaystyle\frac{\partial E^{\alpha_{k}}}{\partial\varphi^{k}}$
$\displaystyle=$ $\displaystyle J_{g}^{\alpha_{k}}\big{(}\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}:\frac{\partial\mathbf{F}_{e}^{\alpha_{k}}}{\partial\varphi^{k}}\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle J_{g}^{\alpha_{k}}\big{(}\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}:\frac{\partial\mathbf{F}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{G}^{\alpha_{k}}]^{-1}\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle J_{g}^{\alpha_{k}}\big{(}\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{M}^{\alpha_{k}}]^{-1}[\mathbf{G}^{\alpha_{k}}]^{-1}\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle J_{g}^{\alpha_{k}}\big{(}\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{M}^{\alpha_{k}}]^{-1}[\mathbf{G}^{\alpha_{k}}]^{-1}\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle J_{g}^{\alpha_{k}}\big{(}\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}[\mathbf{F}^{\alpha_{k}}][\mathbf{G}^{\alpha_{k}}]^{-1}\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle J_{g}^{\alpha_{k}}\big{(}\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}[\mathbf{F}_{e}^{\alpha_{k}}]\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle
J_{g}^{\alpha_{k}}\big{(}\underbrace{\frac{\partial
W^{\alpha_{k}}}{\partial\mathbf{F}_{e}^{\alpha_{k}}}[\mathbf{F}_{e}^{\alpha_{k}}]^{T}}_{J_{e}^{\alpha_{k}}\bm{\sigma}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}\big{)}A_{0}^{\alpha_{k}}$
$\displaystyle=$ $\displaystyle
J_{g}^{\alpha_{k}}{J_{e}^{\alpha_{k}}A_{0}^{\alpha_{k}}\bm{\sigma}^{\alpha_{k}}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}$
$\displaystyle=$ $\displaystyle
A^{\alpha_{k}}\bm{\sigma}^{\alpha_{k}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\det{\mathbf{m}^{\alpha_{k}}}\bm{\sigma}^{\alpha_{k}}:\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\bm{\sigma}^{\alpha_{k}}:\big{(}\det{\mathbf{m}^{\alpha_{k}}}\frac{\partial\mathbf{m}^{\alpha_{k}}}{\partial\varphi^{k}}[\mathbf{m}^{\alpha_{k}}]^{-1}\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\bm{\sigma}^{\alpha_{k}}\varphi^{k\perp,\alpha_{k}}.$

###  6.2 The derivation of Eq.(15) and Eq.(16)

Here, we write $\nabla E(\bm{\varphi}_{i}):=\partial
E/\partial\bm{\varphi}_{i}$ and
$\nabla^{2}E(\bm{\varphi}_{i}):=\partial^{2}E/(\partial\bm{\varphi}_{i}\partial\bm{\varphi}_{i})$
for detailed derivation. First we derive $\Delta s_{stable}$ from Eq.(15) for
the iteration

$\displaystyle\bm{\varphi}_{i+1}-\bm{\varphi}_{i}=-{\Delta s}\nabla
E(\bm{\varphi}_{i}).$ (18)

Assuming $||\bm{\varphi}_{i+1}-\bm{\varphi}_{i}||_{2}\ll 1$, we see from the
Taylor expansion

$E(\bm{\varphi}_{i+1})-E(\bm{\varphi}_{i})=\nabla
E(\bm{\varphi}_{i})(\bm{\varphi}_{i+1}-\bm{\varphi}_{i})+\frac{1}{2}(\bm{\varphi}_{i+1}-\bm{\varphi}_{i})^{T}\nabla^{2}E(\bm{\varphi}_{i})(\bm{\varphi}_{i+1}-\bm{\varphi}_{i})+o(||\bm{\varphi}_{i+1}-\bm{\varphi}_{i}||_{2}^{2})$

that

$E(\bm{\varphi}_{i+1})-E(\bm{\varphi}_{i})\approx-{\Delta s}||\nabla
E(\bm{\varphi}_{i})||_{2}^{2}+\frac{1}{2}{\Delta s}^{2}\nabla
E(\bm{\varphi}_{i})^{T}\nabla^{2}E(\bm{\varphi}_{i})\nabla
E(\bm{\varphi}_{i}).$

To ensure stability $E(\bm{\varphi}_{i+1})-E(\bm{\varphi}_{i})\leq 0$, we have

$\displaystyle\nabla E(\bm{\varphi}_{i})^{T}\big{(}\mathbf{I}-\frac{\Delta
s}{2}\nabla^{2}E(\bm{\varphi}_{i})\big{)}\nabla E(\bm{\varphi}_{i})\geq 0.$
(19)

When $\bm{\varphi}_{i}$ is a local maximizer such that
$\nabla^{2}E(\bm{\varphi}_{i})\leq 0$, Eq.(19) always holds. Otherwise $\Delta
s$ needs to satisfy

$\displaystyle\Delta s\leq{\Delta
s}_{stable}:={2}/{\rho^{+}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}}$

where $\rho^{+}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}$ is the maximum of
the positive eigenvalues of $\nabla^{2}E(\bm{\varphi}_{i})$.

Next we derive $\Delta s_{opt}$ from Eq.(16). First we rewrite the iteration
Eq.(18) as

$\displaystyle\bm{\varphi}_{i+1}-\bm{\varphi}_{*}$ $\displaystyle=$
$\displaystyle\bm{\varphi}_{i}-\bm{\varphi}_{*}-\Delta s\big{(}\nabla
E(\bm{\varphi}_{i})-\nabla E(\bm{\varphi}_{*})\big{)}$ $\displaystyle=$
$\displaystyle\bigg{(}\mathbf{I}-\Delta
s\int_{0}^{1}\nabla^{2}E\big{(}\bm{\varphi}_{i}+r(\bm{\varphi}_{*}-\bm{\varphi}_{i})\big{)}dr\bigg{)}(\bm{\varphi}_{i}-\bm{\varphi}_{*})$
$\displaystyle=$ $\displaystyle\bigg{(}\mathbf{I}-\Delta
s\int_{0}^{1}\nabla^{2}E(\bm{\varphi}_{i})dr\bigg{)}(\bm{\varphi}_{i}-\bm{\varphi}_{*})$
$\displaystyle-\bigg{(}\Delta
s\int_{0}^{1}\big{(}\nabla^{2}E\big{(}\bm{\varphi}_{i}+r(\bm{\varphi}_{*}-\bm{\varphi}_{i})\big{)}-\nabla^{2}E(\bm{\varphi}_{i})\big{)}dr\bigg{)}(\bm{\varphi}_{i}-\bm{\varphi}_{*}).$

where $\bm{\varphi}_{*}$ is a critical point such that $\nabla
E(\bm{\varphi}_{*})=0$. So far we have only assumed that $\nabla
E(\bm{\varphi}_{*})=0$. Following this, we have

$\displaystyle||\bm{\varphi}_{i+1}-\bm{\varphi}_{*}||_{2}$ $\displaystyle\leq$
$\displaystyle||\mathbf{I}-\Delta
s\int_{0}^{1}\nabla^{2}E(\bm{\varphi}_{i})dr||_{2}||\bm{\varphi}_{i}-\bm{\varphi}_{*}||_{2}$
$\displaystyle+||\Delta
s\int_{0}^{1}\big{(}\nabla^{2}E\big{(}\bm{\varphi}_{i}+r(\bm{\varphi}_{*}-\bm{\varphi}_{i})\big{)}-\nabla^{2}E(\bm{\varphi}_{i})\big{)}dr||_{2}||\bm{\varphi}_{i}-\bm{\varphi}_{*}||_{2}.$

Further assuming $\nabla^{2}E$ is Lipschitz continuous with constant $L$, we
have

$\displaystyle||\bm{\varphi}_{i+1}-\bm{\varphi}_{*}||_{2}$ $\displaystyle\leq$
$\displaystyle||\mathbf{I}-\Delta
s\int_{0}^{1}\nabla^{2}E(\bm{\varphi}_{i})dr||_{2}||\bm{\varphi}_{i}-\bm{\varphi}_{*}||_{2}$
$\displaystyle+\frac{\Delta
sL}{2}||\bm{\varphi}_{i}-\bm{\varphi}_{*}||^{2}_{2}.$

Thus, the rate of convergence
$\lambda:=||\bm{\varphi}_{i+1}-\bm{\varphi}_{*}||_{2}/||\bm{\varphi}_{i}-\bm{\varphi}_{*}||_{2}$
is bounded above by

$\displaystyle\lambda$ $\displaystyle\leq$
$\displaystyle\rho\big{(}\mathbf{I}-\Delta
s\nabla^{2}E\big{(}\bm{\varphi}_{i}\big{)}\big{)}+\frac{\Delta
sL}{2}||\bm{\varphi}_{i}-\bm{\varphi}_{*}||_{2}$ (20)

where $\rho(\Box)$ is the spectral radius of $\Box$.

Finally we assume $||\bm{\varphi}_{i}-\bm{\varphi}_{*}||_{2}\ll 1$, the second
term in Eq.(20) can be dropped:

$\displaystyle\lambda\leq\lambda_{bound}(\Delta s)$ $\displaystyle=:$
$\displaystyle\rho\big{(}\mathbf{I}-\Delta
s\nabla^{2}E\big{(}\bm{\varphi}_{i}\big{)}\big{)}$ $\displaystyle=$
$\displaystyle\max{\bigg{(}|1-\Delta
s\lambda_{min}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}|,|1-\Delta
s\lambda_{max}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}|\bigg{)}},$

and our aim is to find $\Delta s_{opt}=\arg\min_{\Delta
s}\lambda_{bound}(\Delta s)\geq 0$.

Notice if $\nabla^{2}E(\bm{\varphi}_{i})$ is not positive semi-definite
($\lambda_{min}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}<0$), we have
$\lambda_{bound}(\Delta s)>1$ for any $\Delta s>0$. Thus we consider
$\nabla^{2}E(\bm{\varphi}_{i})\geq 0$ and find

$\Delta
s_{opt}={2}/\bigg{(}\lambda_{max}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}+\lambda_{min}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}\bigg{)}$

as the minimizer of $\lambda_{bound}(\Delta s)$ (Eq.(16) in the main text),
with the corresponding optimal rate

$\lambda_{opt}=\frac{\lambda_{max}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}-\lambda_{min}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}}{\lambda_{max}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}+\lambda_{min}\big{(}\nabla^{2}E(\bm{\varphi}_{i})\big{)}}.$

Notice the ${\Delta s}_{opt}$ is smaller than ${\Delta s}_{stable}$ from Eq.
(15) when $\nabla^{2}E(\bm{\varphi}_{i})>0$. Thus we choose
$\Delta{s}=\min{(\Delta{s}_{opt},\Delta s_{stable})}$.

## Acknowledgment

This work is partially supported by NSF-DMS-2012330 and CAREER-2144372. The
author expresses gratitude to Gilles Francfort for providing helpful feedback
on the manuscript.

## References

  * Ambrosi et al. [2019] Ambrosi, D., Ben Amar, M., Cyron, C.J., DeSimone, A., Goriely, A., Humphrey, J.D., Kuhl, E., 2019.  Growth and remodelling of living tissues: Perspectives, challenges and opportunities.  Journal of the Royal Society Interface 16\.  doi:10.1098/rsif.2019.0233.
  * Audoly and Pomeau [2000] Audoly, B., Pomeau, Y., 2000\.  Elasticity and geometry, in: Peyresq Lectures on Nonlinear Phenomena. World Scientific, pp. 1–35.
  * Balbi et al. [2015] Balbi, V., Kuhl, E., Ciarletta, P., 2015.  Morphoelastic control of gastro-intestinal organogenesis: Theoretical predictions and numerical insights.  Journal of the Mechanics and Physics of Solids 78, 493–510.
  * Ben Amar [2023] Ben Amar, M., 2023.  Creases and cusps in growing soft matter.  arXiv preprint arXiv:2309.11412 .
  * Ben Amar and Goriely [2005] Ben Amar, M., Goriely, A., 2005\.  Growth and instability in elastic tissues.  Journal of the Mechanics and Physics of Solids 53, 2284–2319.
  * Budday et al. [2015] Budday, S., Steinmann, P., Goriely, A., Kuhl, E., 2015\.  Size and curvature regulate pattern selection in the mammalian brain.  Extreme Mechanics Letters 4, 193–198.
  * Ciarletta and Ben Amar [2012] Ciarletta, P., Ben Amar, M., 2012\.  Growth instabilities and folding in tubular organs: a variational method in non-linear elasticity.  International Journal of Non-Linear Mechanics 47, 248–257.
  * Goriely [2017] Goriely, A., 2017.  The mathematics and mechanics of biological growth.  Springer.
  * Goriely and Ben Amar [2005] Goriely, A., Ben Amar, M., 2005\.  Differential growth and instability in elastic shells.  Physical review letters 94, 198103\.
  * Holzapfel [2002] Holzapfel, G.A., 2002.  Nonlinear solid mechanics: a continuum approach for engineering science.
  * Javili et al. [2015] Javili, A., Dortdivanlioglu, B., Kuhl, E., Linder, C., 2015\.  Computational aspects of growth-induced instabilities through eigenvalue analysis.  Computational Mechanics 56, 405–420.
  * Jin et al. [2011] Jin, L., Cai, S., Suo, Z., 2011.  Creases in soft tissues generated by growth.  Europhysics Letters 95, 64002\.
  * Kröner [1959] Kröner, E., 1959.  Allgemeine kontinuumstheorie der versetzungen und eigenspannungen.  Archive for Rational Mechanics and Analysis 4, 273–334.
  * Lee [1969] Lee, E.H., 1969.  Elastic-plastic deformation at finite strains.  Trans.ASMEJ.Appl.Mech. 54, 1–6.
  * Li et al. [2011a] Li, B., Cao, Y.P., Feng, X.Q., Gao, H., 2011a.  Surface wrinkling of mucosa induced by volumetric growth: theory, simulation and experiment.  Journal of the Mechanics and Physics of Solids 59, 758–774.
  * Li et al. [2011b] Li, B., Jia, F., Cao, Y.P., Feng, X.Q., Gao, H., 2011b.  Surface wrinkling patterns on a core-shell soft sphere.  Physical review letters 106, 234301\.
  * Moulton and Goriely [2011] Moulton, D., Goriely, A., 2011\.  Circumferential buckling instability of a growing cylindrical tube.  Journal of the Mechanics and Physics of Solids 59, 525–537.
  * Papastavrou et al. [2013] Papastavrou, A., Steinmann, P., Kuhl, E., 2013.  On the mechanics of continua with boundary energies and growing surfaces.  Journal of the Mechanics and Physics of Solids 61, 1446–1463.
  * Rodriguez et al. [1994] Rodriguez, E.K., Hoger, A., McCulloch, A.D., 1994.  Stress-dependent finite growth in soft elastic tissues.  Journal of Biomechanics 27, 455–467.  doi:10.1016/0021-9290(94)90021-3.
  * Tallinen et al. [2013] Tallinen, T., Biggins, J.S., Mahadevan, L., 2013.  Surface sulci in squeezed soft solids.  Physical review letters 110, 024302\.
  * Tallinen et al. [2016] Tallinen, T., Chung, J.Y., Rousseau, F., Girard, N., Lefèvre, J., Mahadevan, L., 2016\.  On the growth and form of cortical convolutions.  Nature Physics 12, 588–593.
  * Wu and Ben Amar [2015] Wu, M., Ben Amar, M., 2015.  Modelling fibers in growing disks of soft tissues.  Mathematics and Mechanics of Solids 20, 663–679.