EPJ Web of Conferences Lattice2017

11institutetext: Institut fÃr Theoretische Physik, Universität MÃnster,
Wilhelm-Klemm-Str. 9, D-48149 MÃnster 22institutetext: Theoretisch-
Physikalisches Institut, Universität Jena, Max-Wien-Platz 1, D-07743 Jena
33institutetext: Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, D-22603
Hamburg 44institutetext: Fakultät fÃr Physik, Universität Regensburg,
Universitätsstr. 31, D-93053 Regensburg

# Supermultiplets in $\mathcal{N}$=$1$ SUSY SU(2) Yang-Mills Theory

Sajid Ali 11 Georg Bergner 22 Henning Gerber Speaker, E-mail: h.gerber@uni-
muenster.de11 Pietro Giudice 11 Simon Kuberski 11 István Montvay 33 Gernot
MÃnster 11 Stefano Piemonte 44 Philipp Scior 11

###### Abstract

We study $\mathcal{N}=1$ supersymmetric Yang-Mills theory (SYM) on the
lattice. The non-perturbative nature of supersymmetric field theories is still
largely unknown. Similarly to QCD, SYM is confining and contains strongly
bound states. Applying the variational method together with different smearing
techniques we extract masses of the lightest bound states such as gluino-glue,
glueball and mesonic states. As these states should form supermultiplets, this
study allows to check whether SYM remains supersymmetric also on the quantum
level.

##  1 Introduction

The Standard Model of particle physics accurately describes current particle
physics experiments. However, it is neither theoretically plausible as a
fundamental theory (e.g. hierarchy problem), nor does it explain the
observation of dark matter. It should therefore be viewed as a low energy
effective theory of some more complete theory of physics. Deviations from the
Standard Model are expected to be observed when measuring with more precision
or at higher energies.

One possible extension of the Standard Model which cures the hierarchy problem
and also provides candidates for dark matter is supersymmetry. Supersymmetry
relates bosonic and fermionic fields to one another and implies that for every
particle there exists a supersymmetric partner particle of the same mass and a
spin different by $1/2$. Since this is not observed by current experiments,
supersymmetry must be broken if it exists.

Supersymmetric models have been analyzed extensively within the framework of
perturbation theory. On the other hand, it is desirable to gain more insight
in their non-perturbative properties. We study $\mathcal{N}=1$ supersymmetric
Yang-Mills theory (SYM) on the lattice. It is the simplest supersymmetric
theory containing non-abelian gauge fields. Important topics to be studied
are, besides others, the question whether a supersymmetric continuum limit
exists and what the particle spectrum of such a theory looks like.

When treated within the framework of lattice field theory, supersymmetry is
broken due to the breaking of Poincaré invariance by the lattice. In Ref.
Curci:1986sm  it has been suggested that $\mathcal{N}=1$ SYM can be studied on
the lattice by extrapolating to the chiral and continuum limit, see
Bergner:2016sbv  for a detailed review. An important check of supersymmetry
restoration in this limit is to study the multiplet formation which was
predicted in Ref. Veneziano:1982ah  and Ref. Farrar:1997fn . The authors argue
that if SYM is unbroken, chiral supermultiplets are formed from glueball
states with quantum numbers $0^{++}$ and $0^{-+}$ (in $\text{J}^{PC}$
notation), a-$f_{0}$ and a-$\eta^{\prime}$ mesons111The names stem from the
similarity of these states with the corresponding particles in QCD. The
fermionic constituents are, however, in the adjoint and not in the fundamental
representation of the group, see section 1.1. and additional gluino-glue
states.

The groundstates of these particles have been investigated in Ref.
Bergner:2015adz  for SYM with the gauge group $\mathrm{SU}(2)$. Indeed, a mass
degeneracy of these states (except for the $0^{-+}$-glueball) was found in the
chiral and continuum limit. The goal of this work is to analyze also the
masses of the corresponding first excited states.

###  1.1 The Theory

The Lagrangian of SYM in the continuum reads

$\mathcal{L}=-\frac{1}{4}F_{\mu\nu}^{a}F^{\mu\nu,a}+\frac{i}{2}\bar{\lambda}^{a}\left(\not{D}\lambda\right)^{a}-\frac{m_{\text{g}}}{2}\bar{\lambda}^{a}\lambda^{a},$
(1)

where $F^{\mu\nu,a}$ is the field strength tensor built from the gluon field
$A_{\mu}^{a}$. The corresponding superpartner is the gluino field
$\lambda^{\alpha,a}$. Supersymmetry requires the gluino to be in the adjoint
representation of the gauge group and to be of Majorana type, i.e. fulfill the
Majorana condition

$\bar{\lambda}=\lambda^{T}C,$ (2)

where $C$ is the charge conjugation matrix. The covariant derivative for the
gauge group $\mathrm{SU}(2)$ is given by

$\left(D_{\mu}\lambda\right)^{a}=\partial_{\mu}\lambda^{a}+g\epsilon_{abc}A^{b}_{\mu}\lambda^{c}.$
(3)

The gluino mass term breaks supersymmetry softly. In order to approach the
chiral limit, the bare gluino mass $m_{\text{g}}$ needs to be fine-tuned such
that the renormalized gluino mass vanishes, see also thiscontrib102 .

For the lattice discretization we use Wilson fermions and a tree level
Symanzik improved gauge action. Additionally, we use one level of stout
smearing Morningstar:2003gk  for the links in the fermionic part of the action
in order to improve the signal-to-noise ratio of the measurements.

##  2 Techniques

###  2.1 Variational Method

For the extraction of the excited state masses we use the variational method
as proposed in Ref. Luscher:1990ck . It requires solving a generalized
eigenvalue problem (GEVP) of the form

$C(t)\vec{v}^{\ (n)}=\omega^{(n)}(t,t_{0})C(t_{0})\vec{v}^{\ (n)}.$ (4)

Here $C(t)$ is the correlation matrix, $\omega^{(n)}(t,t_{0})$ and $\vec{v}^{\
(n)}$ are the corresponding generalized eigenvalues and eigenvectors. The
entries of the correlation matrix are given by correlation functions of
different interpolating fields $O_{i}$ of the respective state:

$C_{ij}(t)=\left\langle O_{i}(t)O_{j}^{\dagger}(0)\right\rangle.$ (5)

The large-$t$ behavior of the generalized eigenvalues allows to extract the
mass $m_{n}$ of the excited states by fitting the eigenvalues
$\omega^{(n)}(t,t_{0})$ to exponentials:

$\underset{t\rightarrow\infty}{\lim}\
\omega^{(n)}(t,t_{0})\propto\mathrm{e}^{-m_{n}\left(t-t_{0}\right)}\left(1+\mathcal{O}\left(\mathrm{e}^{-\Delta
m_{n}(t-t_{0})}\right)\right),\quad\Delta m_{n}=\underset{l\neq
n}{\min}\left|m_{l}-m_{n}\right|.$ (6)

For a good signal it is preferable to choose a set of interpolating fields
which on the one hand have good overlaps with the physical states and on the
other hand have only small mutual overlaps.

###  2.2 Interpolating fields

The interpolating fields $O_{i}$ used in the correlation matrix (5) have to
match the quantum numbers of the respective state to be analyzed. The basic
interpolating fields that we use are described below. In order to improve the
signal-to-noise ratio and to create different interpolating fields for each
state we apply different smearing methods to these operators (see Sect. 2.3).

The basic interpolating fields for glueballs are built from gauge link loops
that respect the parity quantum number of the respective state. For the
$0^{++}$-glueballs we use a sum of gauge plaquettes

$O_{\text{gb}^{++}}(x)=\text{Tr}\left[P_{12}(x)+P_{23}(x)+P_{31}(x)\right],$
(7)

where $P_{ij}$ denotes a plaquette in the $i$-$j$ direction (see Fig. 1). For
the pseudo-scalar glueball $0^{-+}$ we use

$O_{\text{gb}^{-+}}(x)=\sum_{R\in\mathbf{O}_{h}}[\text{Tr}\left[\mathcal{C}(x)\right]\text{-}\text{Tr}\left[P\mathcal{C}(x)\right],$
(8)

where the sum is over all rotations of the cubic group $\mathbf{O}_{h}$ and
$P\mathcal{C}$ is the parity conjugate of loop $\mathcal{C}$; it is depicted
in Fig. 1.

Figure 1: Plaquette $P_{ij}(x)$ and gauge-loop $\mbox{$\mathcal{C}$}_{ijk}(x)$
used in the interpolating field of $0^{-+}$-glueball.

The basic interpolating fields for the mesons are

$\displaystyle O_{\text{a-f}_{0}}(x)=$
$\displaystyle\bar{\lambda}(x)\lambda(x),\qquad
O_{\text{a-}\eta^{\prime}}(x)=\bar{\lambda}(x)\gamma_{5}\lambda(x).$ (9)

When inserted into the correlation matrix, Wick contractions of these fields
lead to connected and disconnected pieces

$\langle\bar{\lambda}(x)\Gamma\lambda(x)\bar{\lambda}(y)\Gamma\lambda(y)\rangle=\text{Tr}\left[\Gamma\Delta(x,x)\right]\text{Tr}\left[\Gamma\Delta(y,y)\right]-2\text{Tr}\left[\Gamma\Delta(x,y)\Gamma\Delta(y,x)\right],$
(10)

where $\Delta(x,y)$ denotes the propagator from $x$ to $y$ (spin and group
indices suppressed) and $\Gamma$ represents $\mathbb{1}$ or $\gamma_{5}$.

The interpolating field of the gluino-glue state is given by

$O_{\text{g-g}}^{\alpha}(x)=\sum_{i<j=1}^{3}\sigma_{ij}^{\alpha\beta}\text{Tr}\left[P_{ij}(x)\lambda^{\beta}(x)\right],$
(11)

with $\sigma_{\mu\nu}=\frac{\text{i}}{2}[\gamma_{\mu},\gamma_{\nu}].$

###  2.3 Smearing Methods

In order to improve the signal-to-noise ratio and to obtain different
interpolating fields $O_{i}$ for the variational method, we use APE-smearing
Albanese:1987ds  for the gauge links $U_{\mu}(x)$ and gauge invariant Jacobi
smearing Allton:1993wc ; Best:1997qp  for the gluino field $\lambda(x)$.

An APE-smearing step is given by

$U^{\prime}_{\mu}(x)=U_{\mu}(x)+\epsilon_{\text{APE}}\sum_{\nu=\pm
1,\nu\neq\mu}^{\pm
3}U_{\nu}^{\dagger}(x+\hat{\mu})U_{\mu}(x+\hat{\nu})U_{\nu}(x),$ (12)

where $\epsilon_{\text{APE}}$ is a parameter of the smearing method that needs
to be optimized to maximize the signal-to-noise ratio. The APE-smearing step
can be iteratively applied to obtain higher smearing levels. The gluino field
is smeared by applying the Jacobi smearing function:

$\lambda^{\prime\beta b}(\vec{x},t)=\sum_{\vec{y}}F_{\beta b,\alpha
a}(\vec{x},\vec{y},t)\lambda^{\alpha a}(\vec{y},t),$ (13)

where the smearing function $F_{\beta b,\alpha a}(\vec{x},\vec{y},t)$ is given
by

$F_{\beta b,\alpha
a}(\vec{x},\vec{y},t)=\delta_{\vec{x},\vec{y}}\,\delta_{\beta\alpha}+\delta_{\beta\alpha}\sum_{i=1}^{N_{\text{J}}}\left(\kappa_{\text{J}}\sum_{\mu=1}^{3}\left[\delta_{\vec{y},\vec{x}+\hat{\mu}}U_{\mu}(\vec{x},t)+\delta_{\vec{y}+\hat{\mu}}U_{\mu}^{\dagger}(\vec{x},t)\right]\right)^{i}_{ba},$
(14)

with the smearing parameter $\kappa_{\text{J}}$ and the smearing level
$N_{\text{J}}$.

###  2.4 Identifying the eigenvalues

The errors of the generalized eigenvalues that are used to fit the masses and
the final errors of the masses are estimated by resampling the data using the
Jackknife method.

Fitting an exponential to the $n$-th generalized eigenvalue
$\omega^{(n)}(t,t_{0})$ requires identifying it out of all the generalized
eigenvalues at $t$ on each Jackknife sample. Originally it was proposed to
identify the eigenvalues at each $t$ by their magnitude,
$\omega^{(0)}>\omega^{(1)}>\omega^{(2)},\ldots$ . However, due to limited
statistics the eigenvalues fluctuate on the Jackknife samples. If the
fluctuations of two or more eigenvalues overlap, it is not suitable any longer
to identify them by their magnitude. In this situation we use a different
approach, which is to identify the generalized eigenvalues by their
corresponding generalized eigenvectors $\vec{v}^{\ (n)}$ in the following way:

  1. 1.

Solve the GEVP for all $t$, find a $t_{1}$ where the fluctuations of the
$\omega^{(n)}(t_{1},t_{0})$ are small enough such that they do not overlap.
$t_{1}$ should be chosen as large as possible, so that the eigenvectors have
sufficiently stabilized. Sort the eigenvalues by their magnitude as originally
proposed. The corresponding eigenvectors are called $\vec{w}^{(n)}$.

  2. 2.

Solve the GEVP at all other $t$, call the eigenvectors $\vec{v}^{\
(n)}(t,t_{0})$.

  3. 3.

For each $t$ calculate a matrix of scalar products
$M_{mn}=|\vec{w}^{(m)}\cdot\vec{v}^{\ (n)}(t,t_{0})|$.

  4. 4.

    1. (a)

Scan for largest entry $M_{n_{\text{max}}m_{\text{max}}}$ in $M$.

    2. (b)

Identify $\omega^{(m_{\text{max}})}(t,t_{0})$ to be the eigenvalue
corresponding to $\vec{v}^{(n_{\text{max}})}(t,t_{0})$.

    3. (c)

Delete $m_{\text{max}}$-th row and $n_{\text{max}}$-th column from $M$.

    4. (d)

Repeat 3. until all eigenvalues are identified.

A comparison of this method with the earlier one is shown in Fig. 2.

$0.001$$0.01$$0.1$$1$$0$$2$$4$$6$$8$$10$$12$$14$$t/a$$\omega^{(1)}$$\omega^{(2)}$

$0.0001$$0.001$$0.01$$0.1$$1$$0$$2$$4$$6$$8$$10$$12$$14$$t/a$$\omega^{(1)}$$\omega^{(2)}$

Figure 2: Left: Generalized eigenvalues in the $0^{++}$-channel at
$\beta=1.9$, $\kappa=0.1433$ sorted by magnitude, right: sorted by scalar
products of eigenvectors. Six interpolating fields were used in the GEVP; for
better visibility only two of the six eigenvalues are displayed. Note that the
eigenvalues at $t=5$ have been misidentified when sorted by magnitude. The
identification using scalar products allows to use the data up to $t=8$.

###  2.5 Optimization of the smearing parameters

In order to determine the optimal values for the Jacobi smearing parameters
$\kappa_{\text{J}}$ and $N_{\text{J}}$, we measured the smearing radius
defined as

$R_{\text{J}}^{2}=\dfrac{\sum_{\vec{x}}\left|\vec{x}\right|^{2}\left|F(\vec{x},0)\right|^{2}}{\sum_{\vec{x}}\left|F(\vec{x},0)\right|^{2}}.$

The measurement shows that there is a critical parameter
$0.15<\kappa_{\text{J}}^{\text{c}}<0.2$ (see Fig. 3). For values above
$\kappa_{\text{J}}^{\text{c}}$ the smearing function $(\ref{eq: Jacobi
smearing function})$ diverges in the limit $N_{\text{J}}\rightarrow\infty$,
while for smaller values it converges. Since we use Jacobi smearing to create
new interpolating fields from the basic ones, it is more suitable to choose a
value for $\kappa$ in the diverging regime. In order to keep numerical errors
small, we chose $\kappa_{\text{J}}^{\text{c}}=0.2$ which is just above the
critical $\kappa_{\text{J}}^{\text{c}}$. Optimizing the signal-to-noise ratio,
we chose to use the smearing levels 0, 40 and 80. In principle more smearing
levels could be used, but from our experience it makes the identification of
the eigenvalues (see Sect. 2.4) unstable and therefore does not improve the
results.

A similar analysis was also done for the APE-smearing. In the gluino-glue
measurements we chose $\epsilon_{\text{APE}}=0.4$ and the smearing levels
$0,5,15,25,..,95$. In the glueball measurements we chose the smearing levels
10, 50 and 80.

$0$$1$$2$$3$$4$$5$$6$$7$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$

Radius [a]

Smearing
Level$\kappa_{J}=0.1$$\kappa_{J}=0.15$$\kappa_{J}=0.2$$\kappa_{J}=0.225$$\kappa_{J}=0.25$

Figure 3: Measurement of the Jacobi smearing radius for different values of
$\kappa_{\text{J}}$ at $\beta=1.9,$ and $\kappa=0.14415$. For
$\kappa_{\text{J}}\geq 0.2$ the smearing function diverges for
$N_{\text{J}}\rightarrow\infty$. For $\kappa_{\text{J}}\leq 0.15$ it
converges.

###  2.6 Exploiting mixing

Operators with the same quantum numbers are expected to mix. We expect the
$\text{a-f}_{0}$ meson to mix with the $0^{++}$-glueball since they share the
same quantum numbers. The same is true for the a-$\eta^{\prime}$ meson and the
$0^{-+}$-glueball. This mixing allows both, mesonic and gluonic, interpolators
to be used in the correlation matrix (5) for the extraction of a mixed state.
The full correlation matrix then has the following form:

$C(t)=\left(\begin{array}[]{cc}\left\langle
O_{\text{gb}^{++}}(t)O^{\dagger}_{\text{gb}^{++}}(0)\right\rangle&\left\langle
O_{\text{gb}^{++}}(t)O^{\dagger}_{a-f_{0}}(0)\right\rangle\\\ \left\langle
O_{a-f_{0}}(t)O^{\dagger}_{\text{gb}^{++}}(0)\right\rangle&\left\langle
O_{a-f_{0}}(t)O^{\dagger}_{a-f_{0}}(0)\right\rangle\end{array}\right).$ (15)

Here each entry is a submatrix consisting of the different interpolating
fields for each operator. Since the mesonic and the gluonic operators are very
different from one another, the mutual overlap is expected to be small and
using this larger correlation matrix should improve the signal drastically.

##  3 Preliminary Results

We have used lattices for $\beta=1.9$ at four different hopping parameters
$\kappa\in\\{0.1433,0.14387,0.14415,0.14435\\}$. For each value of $\kappa$
there are more than $10.000$ lattice gauge configurations available of which
only every $8$th configurations has been measured so far. Our results are
therefore only preliminary.

The results (see Fig. 4) indicate that the variational method in combination
with different smearing levels and different interpolating fields can indeed
be used to extract the masses of the excited states. As expected, we observe a
much cleaner signal in the $0^{++}$-channel when using both, mesonic and
gluonic, interpolators than when using them separately. Curiously, in the
$0^{-+}$-channel the correlation matrix is block-diagonal within errors, i.e.
there seems to be no mixing between the a-$\eta^{\prime}$ meson and the
$0^{-+}$-glueball. The groundstate mass of the $0^{-+}$-glueball seems to be
in the range of the excited states. These two findings might be an indication
that the $0^{-+}$-component of the lowest chiral multiplet only consists of
a-$\eta^{\prime}$ and that the $0^{-+}$-glueball actually belongs to a
different, heavier multiplet as was already proposed in Ref. Bergner:2015adz .

For the final results the statistics will be increased and a detailed analysis
of autocorrelations will be done. A chiral extrapolation to the critical
$\kappa_{\text{c}}$ will be performed. Previous results Bergner:2015adz
indicate that the ensemble at $\beta=1.9$ is already close to the continuum
limit. To provide an estimate of lattice discretization errors, we will also
analyze a coarser lattice spacing at $\beta=1.75$.

$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0$$0.02$$0.04$$0.06$$0.08$$0.1$

$am$

$(am_{\pi})^{2}$$0^{++}$-mixed first exciteda-$\eta^{\prime}$ first
excited$0^{-+}$-glueball groundstategluino-glue first excited

Figure 4: Preliminary results for $\beta=1.9$. The extracted masses are
plotted against the measured squared mass of the unphysical adjoint pion
$m_{\pi}$. The gluino mass scales proportionally to $m_{\pi}^{2}$, thus the
chiral limit is at $m_{\pi}^{2}=0$ Munster:2014cja .

## Acknowledgements

The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS)
for providing computing time for a GCS Large-Scale Project on the GCS share of
the supercomputer JUQUEEN at JÃlich Supercomputing Centre (JSC). GCS is the
alliance of the three national supercomputing centres HLRS (Universität
Stuttgart), JSC (Forschungszentrum JÃlich), and LRZ (Bayerische Akademie der
Wissenschaften), funded by the German Federal Ministry of Education and
Research (BMBF) and the German State Ministries for Research of Baden-
WÃrttemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF). Further
computing time has been provided on the supercomputers JURECA at JSC and on
the compute cluster PALMA and NWZPHI of the University of MÃnster. This work
is supported by the Deutsche Forschungsgemeinschaft (DFG) through the Research
Training Group “GRK 2149: Strong and Weak Interactions - from Hadrons to Dark
Matter”.

## References

  * (1) G. Curci, G. Veneziano, Nucl. Phys. B292, 555 (1987)
  * (2) G. Bergner, S. Catterall, Int. J. Mod. Phys. A31, 1643005 (2016), arXiv:1603.04478
  * (3) G. Veneziano, S. Yankielowicz, Phys. Lett. 113B, 231 (1982)
  * (4) G. R. Farrar, G. Gabadadze, M. Schwetz, Phys. Rev. D58, 015009 (1998), arXiv:hep-th/9711166
  * (5) G. Bergner, P. Giudice, I. Montvay, G. MÃnster, S. Piemonte, JHEP 03, 080 (2016), arXiv:1512.07014
  * (6) S. Ali, G. Bergner, H. Gerber, P. Giudice, I. Montvay, G. MÃnster, S. Piemonte, P. Scior, _Ward identities in N = 1 supersymmetric SU(3) Yang-Mills theory on the lattice_ , in _Proceedings,35th International Symposium on Lattice Field Theory (Lattice2017): Granada, Spain_, to appear in EPJ Web Conf., arXiv:1710.01000
  * (7) C. Morningstar, M. J. Peardon, Phys. Rev. D69, 054501 (2004), arXiv:hep-lat/0311018
  * (8) M. LÃscher, U. Wolff, Nucl. Phys. B339, 222 (1990)
  * (9) M. Albanese et al. (APE Collaboration), Phys. Lett. B192, 163 (1987)
  * (10) C. R. Allton et al. (UKQCD Collaboration), Phys. Rev. D47, 5128 (1993), arXiv:hep-lat/9303009
  * (11) C. Best et al., Phys. Rev. D56, 2743 (1997), arXiv:hep-lat/9703014
  * (12) G. MÃnster, H. StÃwe, JHEP 05, 034 (2014), arXiv:1402.6616


# Buckling without bending morphogenesis:
Nonlinearities, spatial confinement, and branching hierarchies

M. C. Gandikota<EMAIL_ADDRESS>Physics Department and BioInspired
Institute, Syracuse University, Syracuse, NY 13244 USA  J. M. Schwarz
<EMAIL_ADDRESS>Physics Department and BioInspired Institute,
Syracuse University, Syracuse, NY 13244 USA  Indian Creek Farm, Ithaca, NY
14850 USA

(June 29, 2021)

###### Abstract

During morphogenesis, a featureless convex cerebellum develops folds. As it
does so, the cortex thickness is thinnest at the crest (gyri) and thickest at
the trough (sulci) of the folds. This observation cannot be simply explained
by elastic theories of buckling. A recent minimal model explained this
phenomenon by modeling the developing cortex as a growing fluid under the
constraints of radially spanning elastic fibers, a plia membrane and a
nongrowing sub-cortex (Engstrom, et. al., PRX 2019). In this minimal buckling
without bending morphogenesis (BWBM) model, the elastic fibers were assumed to
act linearly with strain. Here, we explore how nonlinear elasticity influences
shape development within BWBM. The nonlinear elasticity generates a quadratic
nonlinearity in the differential equation governing the system’s shape and
leads to sharper troughs and wider crests, which is an identifying
characteristic of cerebellar folds at later stages in development. As
developing organs are typically not in isolation, we also explore the effects
of steric confinement, and observe flattening of the crests. Finally, as a
paradigmatic example, we propose a hierarchical version of BWBM from which a
novel mechanism of branching morphogenesis naturally emerges to qualitatively
predict later stages of the morphology of the developing cerebellum.

##  I Introduction

When morphogenesis is viewed through the lens of a physicist, one of the
typical mechanistic routes to take is morphoelasticity induced by varying
internal stresses Ben Amar  _et al._ (2019). Shape change in response to
internal stresses, within elastic objects which tend to retain shape, provides
us with a large collection of shapes. Differential growth can be a source of
such an internal stress. This reasonable viewpoint is made manifest in the
widespread use of the Euler buckling instability in purely elastic materials
to explain the onset of folds in morphogenesis problems Nelson (2016) as
diverse as cerebra Budday  _et al._ (2015); Richman  _et al._ (1975); Raghavan
_et al._ (1997); Bayly  _et al._ (2013); Tallinen  _et al._ (2014, 2016),
intestinal crypts and villi Hannezo  _et al._ (2011); Shyer  _et al._ (2013),
airway mucus wrinkles Wiggs  _et al._ (1997); Li  _et al._ (2011), tooth
ridges Osborn (2008) and hair-follicle patterns Shyer  _et al._ (2017). Recent
work accounting for cerebellar foliation falls under the same framework
Lejeune  _et al._ (2016, 2019). These models typically predict a
characteristic length scale between folds and state quantitative agreement
between prediction and observation to validate such an approach.

On the other hand, tissue fluidity in early stages of developing embryos has
emerged as a driver of shape change in both zebrafish and the insect Tribolium
castaneum Mongera and et al. (2018); Jain and et al (2020). Fluids, unlike
elastic solids, do not retain their shape. Then how does fluidity - the
antithesis of elasticity - affect the shape of the developing organ?
Presumably biological systems have figured out ways to combine elements of
elasticity and fluidity to bring about an even more complex collection of
shapes at later stages of development. The buckling without bending
morphogenesis (BWBM) model is one such model affirming the cleverness of
biological systems Engstrom  _et al._ (2018). Its origin is rooted in
explaining recent quantitative observations of the developing murine
cerebellum, a much less studied component of the brain as compared to the
cerebrum.

Vertebrate brains are typically divided into three different sectors – the
forebrain, the midbrain, and the hindbrain. While the forebrain consists of
the cerebrum, the hypothalamus, and the thalamus, the hindbrain contains the
spinal cord, medulla oblongata, and the cerebellum. Since the shape of the
forebrain varies more across species than the midbrain and hindbrain, its
development has been the focus of much study. On the other hand, by studying
conserved patterns in the hindbrain, we gain complementary insight into the
morphogenetic processes of the brain.

To be specific, let us characterize the difference in shape between the
cerebrum and cerebellum. First, the shape of the cerebellum has an approximate
cylindrical symmetry, whereas such a symmetry is absent in the cerebrum.
Second, mammalian cerebellums of all sizes produce folds Larsell (1967),
whereas small mammalian cerebrums do not show folds Herculano-Houzel (2009).
Given that the overall size of both organs increases with increasing species
size, it is possible that the differences in shape arises via different
physical shape change mechanisms. If so, do different physical shape change
mechanisms lead to different emergent functionalities of the two? Moreover,
the conservation of the number of the 8-10 primary lobes of the cerebellum
across species demands a treatment on its own footing as it suggests a scale
invariant shape change mechanism. Number of secondary folds however, can
differ across species.

Recent experimental observations of the developing murine cerebellum found
that the proliferating cells in the cerebellar cortex are motile with neighbor
exchanges on the order of minutes Lawton  _et al._ (2019). It was also
observed that the developing cerebellar cortex varies in thickness with the
trend being that the cortex is thinnest at the crest (gyri) and thickest at
the trough (sulci). Moreover, the developing cerebellum is under tension and
not under compression as evidenced by both radial and circumferential cuts.
While these observations cannot be explained by elastic wrinkling theories,
all three of these findings can be explained by the linear BWBM model in which
a growing cerebellar cortex is fluid-like and the sub-cortex is a nongrowing
core Engstrom  _et al._ (2018). The cortex is under tension due to the
presence of radial glial cells spanning the cerebellum as well as the pial
membrane. Finally, Bergmann glial cells spanning the cortex attempt to
maintain constant thickness of the cortex. See Fig. 1. Cell growth in the
presence of such constraints drives a featureless convex cortex to form folds.
In addition, the BWBM model also offers an explanation for the length scale
invariance of the formation of cerebellum folds to understand the conservation
of 8-10 primary lobes across vertebrate species spanning a range of sizes
Engstrom  _et al._ (2018).

The linear BWBM model provides a quantitative framework for the onset of shape
change that manifest as smooth cortex oscillations in the developing
cerebellum. However, as the shape of the cerebellum continues to evolve, we
observe cusped sulci and wide gyri (see Fig. 2). The linear BWBM model, which
can be mapped to a forced, simple harmonic oscillator, cannot account for such
nonlinear phenomena. These observations necessitate the exploration of
nonlinear elasticity within the BWBM model, particularly in the context of
tensioned radial glial cells given that robust nonlinearities have been
observed in stretched cells Fernandez  _et al._ (2006). In addition to
exploring the effect of nonlinear springs in the BWBM model, we will also
explore the impact of spatial confinement on shape change, which imposes a
different form of nonlinearity in the BWBM model. Finally, since we are
exploring stages of shape development beyond the onset of shape change, we
note the hierarchy of folds that emerges in the developing cerebellum. This
hierarchy manifests itself as folds within folds. As the developing cerebellum
grows in size, this leads to a type of branching morphogenesis. While
branching morphogenesis has been thought of in the context of developing
lungs, kidneys, and other organs Varner and Nelson (2014); Hannezo  _et al._
(2017), a hierarchical version of BWBM, given its scale invariant solutions,
naturally emerges as a potential candidate to explain branching morphogenesis
within the cerebellum.

The outline of the manuscript is as follows. In Section II, we review the
linear BWBM model focusing on the early stages of shape change in the
developing cerebellum. Section III discusses the effect of nonlinear springs
that describe the radial fibrous glial cells and its role in sharpening the
sinusoidal sulci of the linear model. Section IV discusses the effect of
steric hindrances on the flattening of gyri and Section V discusses the
hierarchy of subsequent folds emerging at even later stages of cerebellar
development. Section VI summarizes the work and addresses its implications.

##  II Linear buckling without bending morphogenesis model

Figure 1: Schematic of the buckling without bending morphogenesis (BWBM)
model. The radial glial cells span the cerebellum and Bergmann glial cells
span the cortex and resist thickness variations in the cortex. The area of the
sub-cortex does not grow as fast as the cortex and so it is approximated as
nongrowing over some time scale.

Within a timespan of a day, the featureless, convex developing cerebellum
begins to develop folds. The BWBM model provides a physical basis for the
onset of these folds. Unlike the cerebrum, the cylindrical symmetry of the
cerebellum allows for two-dimensional modeling of its sagittal cross section.
The two-dimensional BWBM model Engstrom  _et al._ (2018) offers a quasi-static
description of the foliation i.e., it does not describe the dynamics of growth
but determines the final equilibrium shape for a given set of parameters. As
the parameters change with time, so does the shape. In polar coordinates, the
radius of the cerebellum and the thickness of the cortex, otherwise known as
the external granular layer (EGL), are represented with $r,t$ respectively
(see Fig. 1). The energy functional for the bi-layer model with $\theta$
parametrizing the two degrees of freedom is,

$E\left[r,t,\frac{dt}{d\theta}\right]=\int
d\theta\left\\{k_{r_{0}}(r-r_{0})^{2}-k_{t}(t-t_{0})^{2}+\beta\left(\frac{dt}{d\theta}\right)^{2}\right\\}.$
(1)

Here, all the constants $k_{r_{0}},k_{t},\beta,r_{0},t_{0}$ are positive. The
first term represents the elastic contribution of the radial glial fibrous
cells spanning the cerebellum and the elastic pial membrane surrounding the
cerebellum. Moreover, $k_{r}$ is the elasic modulus for this term with $r_{0}$
as its rest length radius. The second term is a growth potential for the
cortex with $k_{t}$ controlling the contribution of the term and $t_{0}$
setting the rest length thickness. The third term represents the Bergmann
glial fibrous cells resisting thickness variations of the cortex with $\beta$
being the accompanying elastic modulus. Note that the third term is not a
bending energy term. A bending energy term would involve the curvature as
given by the second derivative with $\theta$, unlike the first derivative used
here. Finally, given the slower growth of the sub-cortex/core, as compared to
the cortex, we demand that the area of the sub-cortex does not change, at
least over the time scale of the onset of the foliation, i.e. a day. In other
words,

$\frac{1}{2}\int d\theta(r-t)^{2}=A_{0}=\text{constant}.$ (2)

It is the area conservation that sets up the competition between the radial
elastic energy of the glial fibers and pial membrane with the growth potential
term.

The extremum of the variational problem subject to this constraint yields a
pair of coupled Euler-Lagrange equations for $r,t$. The solution to these
equations are linear sinusoidal functions of the form,

$\begin{split}t(\theta)&=T\;\text{sin}(n\theta+\phi),\\\
r(\theta)&=\left(\frac{1}{1-\epsilon}\right)r_{0}(\theta)-\left(\frac{\epsilon}{1-\epsilon}\right)t(\theta),\end{split}$
(3)

where

$\displaystyle n$ $\displaystyle=\sqrt{\rho(1+\epsilon c/(1-\epsilon))},$
$\displaystyle T$
$\displaystyle=\sqrt{2}\epsilon\sqrt{\frac{A_{0}}{\pi}-\left(\frac{a-t_{0}}{1-\epsilon+\epsilon
c}\right)^{2}}.$

The dimensionless constants are $c:=k_{r_{0}}/k_{t},\epsilon:=\mu/k_{r_{0}}$
and $\rho:=k_{t}/\beta$.

The linear BWBM model successfully explains the cortex thickness oscillations
being out of phase with the sub-cortex deformation for a developing cerebellum
which leads to a thinner cortex at the crest and a thicker cortex at the
trough. This model also allows the amplitude of the cortex thickness
oscillation to be the same size or even greater than the amplitude of the
substrate oscillations. Such phenomena is in contrast to elastic wrinkling
models which fail to wrinkle for thicker cortices at the trough Lagrange  _et
al._ (2016).

Figure 2: Visual comparison of experiment and model. The morphogenesis of
cerebellum involves the development of folds in an initially featureless
convex cerebellum. The top row shows midsagittal cross-sections of mouse
cerebellum prior to birth. Scale bar is 200 $\mu$m. The arrowheads are from
the original figures in Ref. Marzban (2017) and indicate fissures. Reprinted
with permission from Springer Nature. The linear BWBM model reproduces the
smooth sinusoidal like folds that occur at the onset of folding. The BWBM
model with nonharmonic radial glial springs reproduces the sharp cusped sulci
that we see at the later stages of cerebellar development. Note that time is
not explicit in the quasi-static BWBM. We start from system parameters
describing the unstable state and minimize the energy while obeying
constraints to arrive at the nonlinear BWBM configuration. We use brain
folding terminology - gyri (sulci) and oscillator terminology - crests
(troughs) interchangeably to refer to hills (valleys) of oscillations
respectively. Parameters:
$c=0.1,0.155,\epsilon=0.6,\rho=31.3,\tilde{t}_{0}=0.1,0.37,\tilde{k}_{r_{1}}=0,-0.97,\tilde{k}_{r_{2}}=0,1.05$
for the linear and nonlinear BWBM respectively.
[$\tilde{r},d\tilde{r}/d\theta]_{\theta=0}\mbox{ is }[2.01,-0.83]$.

##  III nonlinear elasticity

While the linear BWBM model addresses the onset of shape change, a next follow
up question to ask is – what are the limitations of the linear BWBM model in
explaining the more dramatic shape changes in the developing cerebellum at
later stages of development? As the radial glial cells and the pial membrane
become more stretched by the developing crests, the enhanced stretching may
lead to detachment of the radial glial cells from the pial membrane or may
lead to nonlinear elastic effects. It has long been known that stretched cells
act as nonlinear springs Fernandez  _et al._ (2006) and that collagen, which
is one of dominant fibrous proteins constituting the pial membrane exhibits
nonlinear elasticity as strain increases Storm  _et al._ (2005). Such effects
are not addressed in the linear BWBM model. It would therefore be prudent to
examine the role of nonlinear elasticity in the BWBM model. To do so, we
promote the radial spring constant $k_{r_{0}}$ to $k_{r}(r)$, a radially
dependent spring ‘constant’ -

$k_{r}(r):=k_{r_{0}}+k_{r_{1}}(r-r_{0})+k_{r_{2}}(r-r_{0})^{2}.$ (4)

The above three terms correspond to quadratic, cubic and quartic energy terms
of the radial glial spring potential energy. The cubic energy term brings an
asymmetry in the radial spring energy across $r_{0}$ and the quartic energy
term counteracts the destabilizing effect of the cubic term. The more general
form of the uncoupled differential equation where every spring constant is
allowed to be nonlinear is explored in Appendix VII.1. We nondimensionalize
the problem as
$\tilde{E}:=E/(k_{r_{0}}r_{0}^{2}),\tilde{r}:=r/r_{0},\tilde{t}:=t/r_{0},\tilde{t}_{0}:=t_{0}/r_{0}$
and

$\tilde{K}_{r}(r):=\frac{k_{r}(r)}{k_{r_{0}}}=1+\tilde{k}_{r_{1}}(\tilde{r}-1)+\tilde{k}_{r_{2}}(\tilde{r}-1)^{2},$
(5)

where $\tilde{k}_{r_{1}}=k_{r_{1}}r_{0}/k_{r_{0}}$ and
$\tilde{k}_{r_{2}}=k_{r_{2}}r_{0}^{2}/k_{r_{0}}$. The nondimensional energy
functional with dimensionless variables and coefficients is,

$\begin{split}\tilde{E}\left[\tilde{r},\tilde{t},\frac{d\tilde{t}}{d\theta}\right]=\int
d\theta\bigg{\\{}&\tilde{K}_{r}(r)(\tilde{r}-1)^{2}-\frac{1}{c}(\tilde{t}-\tilde{t}_{0})^{2}\\\
&+\frac{1}{\rho
c}\left(\frac{d\tilde{t}}{d\theta}\right)^{2}\bigg{\\}},\end{split}$ (6)

whose minimization is subject to the constraint,

$\frac{1}{2}\int
d\theta(\tilde{r}-\tilde{t})^{2}=\frac{A_{0}}{r_{0}^{2}}=\text{dimensionless
constant}.$ (7)

The variational problem at hand is,

$\delta\left[\tilde{E}-\epsilon\int d\theta(\tilde{r}-\tilde{t})^{2}\right]=0$
(8)

The corresponding Euler-Lagrange equations are

$\tilde{K}_{r}(r)(\tilde{r}-1)-\epsilon(\tilde{r}-\tilde{t})+\frac{\tilde{K}_{r}^{\prime}(r)}{2}(\tilde{r}-1)^{2}=0,$
(9)

$\frac{1}{\rho
c}\frac{d^{2}\tilde{t}}{d\theta^{2}}+\frac{1}{c}(\tilde{t}-\tilde{t}_{0})-\epsilon(\tilde{r}-\tilde{t})=0.$
(10)

For the purpose of numerically solving these coupled differential equations,
we may as well stop here. However, the exercise of uncoupling the differential
equations points to an important source of nonlinearity. Towards finding it,
we differentiate Eq. 9 twice with $\theta$ to arrive at

$\begin{split}&\left[\tilde{K}_{r}(r)+2(\tilde{r}-1)\tilde{K}_{r}^{\prime}(r)+\frac{(\tilde{r}-1)^{2}}{2}\tilde{K}_{r}^{\prime\prime\prime}(r)-\epsilon\right]\left(\frac{d^{2}\tilde{r}}{d\theta^{2}}\right)\\\
+&\left[3\tilde{K}_{r}^{\prime}(r)+3(\tilde{r}-1)\tilde{K}_{r}^{\prime\prime}(r)+\frac{(\tilde{r}-1)^{2}}{2}\tilde{K}_{r}^{\prime\prime\prime}(r)\right]\left(\frac{d\tilde{r}}{d\theta}\right)^{2}\\\
=&-\epsilon\frac{d^{2}\tilde{t}}{d\theta^{2}}.\end{split}$ (11)

Here we see the nonlinear term $(d\tilde{r}/d\theta)^{2}$. This is a
consequence of the inhomogeneous radial spring constant and the coupling
between the Euler-Lagrange equations brought about by the Lagrange multiplier
in Eq. 8. Solving Eq. 9 for $\tilde{t}$, we have,

$\tilde{t}=\tilde{r}-\frac{1}{\epsilon}\left[(\tilde{r}-1)\tilde{K}_{r}(r)+\frac{\tilde{K}_{r}^{\prime}(r)}{2}(\tilde{r}-1)^{2}\right].$
(12)

Substituting Eq. 12 in Eq. 10 leads to,

$\begin{split}\frac{d^{2}\tilde{t}}{d\theta^{2}}=\rho\bigg{[}&\bigg{(}\bigg{\\{}(\tilde{r}-1)\tilde{K}_{r}(r)+\frac{\tilde{K}_{r}^{\prime}(r)}{2}(\tilde{r}-1)^{2}\bigg{\\}}-\tilde{r}\bigg{)}\\\
\times&\bigg{(}\frac{1+\epsilon c}{\epsilon}\bigg{)}+\epsilon
c\tilde{r}+\tilde{t}_{0}\bigg{]}.\end{split}$ (13)

Substituting, Eq. 13 in Eq. 11, we find the shape equation of the system to
be,

$\begin{split}&~{}~{}~{}\left[\tilde{K}_{r}(r)+2(\tilde{r}-1)\tilde{K}_{r}^{\prime}(r)+\frac{(\tilde{r}-1)^{2}}{2}\tilde{K}_{r}^{\prime\prime}(r)-\epsilon\right]\left(\frac{d^{2}\tilde{r}}{d\theta^{2}}\right)\\\
&+\left[3\tilde{K}_{r}^{\prime}(r)+3(\tilde{r}-1)\tilde{K}_{r}^{\prime\prime}(r)+\frac{(\tilde{r}-1)^{2}}{2}\tilde{K}_{r}^{\prime\prime\prime}(r)\right]\left(\frac{d\tilde{r}}{d\theta}\right)^{2}\\\
&-\rho\epsilon\tilde{r}\\\ &+\rho(1+\epsilon
c)\left[(\tilde{r}-1)\tilde{K}_{r}(r)+\frac{(\tilde{r}-1)^{2}}{2}\tilde{K}_{r}^{\prime}(r)\right]\\\
&+\rho\epsilon t_{0}=0.\end{split}$ (14)

The importance of the term $(d\tilde{r}/d\theta)^{2}$, a source of
nonlinearity in Eq. 14, lies in the robustness of its appearance. Irrespective
of the form of nonlinearity in $k_{r}(r)$, we retain this nonlinear term
whereas the coefficients in Eq. 14 correspondingly vary (see Appendix VII.1).
Such a nonlinearity is also seen when one attempts to uncouple the Lotka-
Volterra equations Davis (1961).

We use the RK45 method of scipy.integrate package in python for numerical
integration of all differential equations in this paper. For generating the
phase-portraits in Fig. 7, we use XPPAUT Ermentrout (2002).

Results of numerical integration in Fig. 3a,b, show asymmetric oscillations -
sharper sulci and smooth gyri. This holds true even with
$\tilde{k}_{r_{1}}=0$, i.e. with no explicit imposed asymmetry in the energy
functional (see Eq. 6). This points to the possible role of the
$(d\tilde{r}/d\theta)^{2}$ nonlinearity in bringing about asymmetric
oscillations. In the context of explaining the sharp sulci in the cerebral
cortex of larger mammals, the sulcification instability in nonlinear elastic
materials has been used Tallinen  _et al._ (2014). Here we observe, sans an
elastic instability, sharper sulci in comparison to their sinusoidal
counterparts. Even under the influence of these nonlinearities, the system
robustly retains the thicker cortical troughs and thinner cortical crests
observed in the linear BWBM model. For small
$\tilde{k}_{r_{1}},\tilde{k}_{r_{2}}$, the number of folds do not change
appreciably in comparison to the linear BWBM model.

Figure 3: Asymmetry in the widths of gyri (crests) and sulci (troughs). a,b)
Numerical solutions in polar and Cartesian coordinates for BWBM with
nonharmonic radial glial springs. Nonsinusoidal sharp sulci can be seen for
systems with nonzero $\tilde{k}_{r_{1}},\tilde{k}_{r_{2}}$ (green, blue,
yellow). c) Phase space orbits of solutions in (a,b). The left-right asymmetry
in the orbits seen as a steeper drop on the left translates to sharper troughs
in the coordinate space. d) Numerical solution in Cartesian coordinates for
the ADO oscillator. Sharper troughs are observed due to the addition of the
nonlinear force term $(dx/dt)^{2}$ to the simple harmonic oscillator. The
acceleration of the oscillator is zero at the points marked with the asterisk.
The width of crest and trough is shown. e) Numerical and analytical solution
of the phase-space orbit for the ADO. The velocity nullcline, in red, is the
locus of points with zero acceleration. The position coordinate (marked in
asterisk) at which the nullcline intersects the orbit decides the widths of
crests and troughs in (d). The shape of the orbit is seen to be especially
similar to the orbits in (c) for the cases of $\tilde{k}_{r_{1}}=0,0.19$. f)
Crest - trough asymmetry for the ADO and the BWBM model as a function of their
respective tuning parameters: $\Gamma$ and
$-\tilde{k}_{r_{1}}=\tilde{k}_{r_{2}}$. The dashed black line is the upper
limit of the asymmetry measure. Parameters 5-lobe BWBM in (a)-(c) and (f):
$c=0.1809,\epsilon=0.9,\rho=15.6,\tilde{t}_{0}=0.7$. For periodicity in (a),
parameter $c$ is appropriately tuned for choices of $\tilde{k}_{r_{1,2}}$.
6-lobe BWBM in (f): $c=0.1,\epsilon=0.6,\rho=31.3,\tilde{t}_{0}=0.5$.
$[\tilde{r},d\tilde{r}/d\theta]_{\theta=0}$ is [1.1,0],[1.25,0] for 5-/6-
lobe.

###  III.1 Assisting-dampening oscillator

In Sec. III, the nonharmonic radial glial springs generate a shape equation 14
which has several sources of nonlinearity. Given the assured presence of the
quadratic nonlinearity $(d\tilde{r}/d\theta)^{2}$ in the shape equation
irrespective of the nonlinearity introduced, we seek to isolate the effect it
has in determining the shape of the system. Towards that end, we study the
simple harmonic oscillator with a velocity dependent force in this section.

The differential equation of motion for a one dimensional, simple harmonic
oscillator of mass $m$ with an external nonlinear forcing term,
$-\gamma(d\tilde{x}/dt)^{2}$ is written as,

$m\frac{d^{2}\tilde{x}}{d\tilde{t}^{2}}=-k\tilde{x}-\gamma\left(\frac{d\tilde{x}}{d\tilde{t}}\right)^{2}.$

The quadratic term assists/dampens the negative/positive velocity respectively
bringing an asymmetry in the problem. We refer to this unphysical oscillator
as the ‘assisting-dampening oscillator’ (ADO). The parameters in this equation
and the length scale $x_{0}$ from the initial condition
$\tilde{x}(0)=\tilde{x}_{0}$ can be used to introduce nondimensional
variables, $x:=\tilde{x}/x_{0}$ and $t:=\tilde{t}/\sqrt{m/k}$. The equation of
motion effectively determined by a single tuning parameter
$\Gamma:=\gamma\tilde{x}_{0}/m$ is,

$\frac{d^{2}x}{dt^{2}}=-x-\Gamma\left(\frac{dx}{dt}\right)^{2}.$ (15)

The nonlinear forcing term is nonconservative since it cannot be written down
as the gradient of a position dependent potential. This problem can also not
be formulated within the framework of Lagrangian mechanics as this velocity
dependent force can not be represented by a function $U(x,\dot{x})$ such that,

$-\Gamma\dot{x}^{2}=-\frac{\partial U}{\partial
x}+\frac{d}{dt}\left(\frac{\partial U}{\partial\dot{x}}\right),$

where $\dot{x}=dx/dt$ Goldstein  _et al._ (2002). However, the time reversal
symmetry of Eq. 15 ensures all trajectories sufficiently close to the
equilibrium point, ($x_{0},y_{0})=(0,0)$ to be closed orbits in the phase
space Strogatz (2018). This translates to periodic motion in the position-time
space (see Fig. 7).

The second order differential equation 15 can be represented as two coupled
first order differential equations,

$\begin{split}\frac{dx}{dt}&=y,\\\ \frac{dy}{dt}&=-x-\Gamma y^{2}.\end{split}$
(16)

Dividing the set of equations, we have,

$y\frac{dy}{dx}=-x-\Gamma y^{2}.$ (17)

In contrast to the differential equation of motion in position-time space (see
Eq. 15), the phase-space orbit is represented by a first-order differential
equation. This makes finding an exact solution to the orbit in phase-space
simpler. The variable transformation, $u=y^{2}$ linearizes the above
differential equation as,

$\frac{du}{dx}+2\Gamma u=-2x.$ (18)

The solution to the homogeneous differential equation,

$\frac{du}{dx}+2\Gamma u=0,$ (19)

is $u(x)=A\;e^{-2\Gamma x}$ where $A$ is a constant that needs to be fixed by
the initial condition. The particular solution to Eq. 18 is
$u(x)=-x/\Gamma+1/(2\Gamma^{2})$. Reverting to $y$, the total solution
satisfying the initial condition $y(-1)=0$ is,

$y(x)=\frac{1}{\sqrt{2}\Gamma}\sqrt{1-2\Gamma
x-(1+2\Gamma)e^{-2\Gamma(x+1)}}.$ (20)

In the limit of $\Gamma\rightarrow 0$, we have $y(x)\sim\pm\sqrt{1-x^{2}}$
which are semi-circular arcs as expected for a simple harmonic oscillator.

The numerical integration of the differential equation is presented in Fig. 3.
The sharp troughs in the position-time space are solely due to the effect of
the quadratic nonlinear term $(dx/dt)^{2}$.

The red lines in Figs. 3e, 7 are velocity nullclines which are the locus of
points in phase-space where the acceleration, $dy/dt=d^{2}x/dt^{2}=0$. From
Eq. 15, it follows that the nullcline is a quadratic function in $x$,

$x=-\Gamma\left(\frac{dx}{dt}\right)^{2}.$ (21)

It is interesting to note that the quadratic nullcline (see Fig. 7b) bears
resemblance to a similarly curved nullcline of the BWBM with nonharmonic
springs (see Fig. 7d). The nullcline’s curvature is dictated by the same
quadratic nonlinear term $(dx/dt)^{2}$ in both cases.

For comparison with this unphysical oscillator, we also study a conservative
Hamiltonian system and observe that it is also capable of exhibiting
oscillations with sharp troughs. The simplest such system is obtained by
adding cubic and quartic potential energy terms to the simple harmonic
oscillator. The former term provides an explicit asymmetry in the energy
functional, and the latter term ensures stability about the equilibrium point.
In Appendix VII.2, this system is explored numerically and analytically using
a perturbation series constructed via the Lindstedt-Poincaré method Jordan
_et al._ (2007).

Figure 4: Space constraint locally flattens the lobes of the linear BWBM
model. a,b) Polar and cartesian representation of the solutions to the
differential equation of the linear BWBM model under a
$\text{tanh}(\tilde{r}-\tilde{r}_{c})$ space constraint (see Eq. 26). The
black confining circles in (a) are the confining walls and are represented by
dashed lines in (b). c) Phase-space orbits, unlike the nonlinear BWBM model
and the ADO model, remain symmetric. Parameters:
$c=0.07,\epsilon=0.9,\rho=15.6,\tilde{K}_{c}=30,q=10^{4},\tilde{t}_{0}=7q,\tilde{r}_{0}=1q$.
[$\tilde{r},\frac{d\tilde{r}}{d\theta}]_{\theta=0}\mbox{ is
}[\tilde{r}_{c},0]$.

###  III.2 A measure for crest - trough asymmetry

Stokes, in his study of propagating waves approaching the shore, discussed the
development of narrow crests and wider troughs Stokes (1880). In this context,
measures for velocity and acceleration skewness have been proposed and
studied. In the folds of cerebellum, which is a stationary spatial
oscillation, a similar asymmetry presents itself in the form of wide
gyri/crest and sharp sulci/trough. We propose the following measure for
quantifying the asymmetry in widths -

$\text{crest - trough
asymmetry}:=\frac{t_{\text{crest}}-t_{\text{trough}}}{t_{\text{crest}}+t_{\text{trough}}}.$
(22)

Here, the length the parameter $\theta$ traverses between consecutive pair of
points at which $d^{2}x/dt^{2}=0$ on the climbing (falling) wave, and the
immediately following falling (climbing) wave defines $t_{\text{crest}}$
($t_{\text{trough}}$), respectively (see Fig. 3d). These points are the points
of intersection of the phase-space orbit with the velocity nullcline (see Fig.
3e). The horizontal mouthed parabolic-shaped nullclines, thus influence the
position of these points of intersection and play a prominent role in bringing
the asymmetry between the widths of the crests and troughs. For the ADO, the
velocity nullcline is exactly parabolic (see Eq. 21 and Fig. 7). The parabolic
shape of the nullcline is due to the quadratic nonlinearity $(dx/dt)^{2}$ in
Eq. 15. For the BWBM model with nonlinear radial glial springs, the nullcline
will not be exactly parabolic due to the presence of other nonlinearities in
Eq. 14.

The measure in Eq. 22 is bounded within $(-1,1)$. For a sinusoidal wave, which
is perfectly symmetric, the measure equals zero. For the case of studying this
asymmetry in the lobes of the nonlinear BWBM model, we use $\theta$ in lieu of
$t$ in Eq. 22. Studying this asymmetry for the ADO and the nonlinear BWBM
model for a range of tuning parameters, we see that the asymmetry scales as
1.04 $\pm$ 0.05 and the nonlinear BWBM model shows a parameter dependent
scaling of 0.61 $\pm$ 0.02 and 0.79 $\pm$ 0.08 for two chosen sets of
parameters of 5 lobed and 6 lobed systems respectively (see Fig. 3f).

##  IV Spatial confinement

A cerebellum does not grow in isolation. It encounters steric effects from the
cerebrum, brain-stem, and the skull. The effect of the skull on a growing
cerebrum has earlier been studied computationally Nie  _et al._ (2010). We,
therefore, are compelled to explore the effects of steric confinement with the
BWBM model, particularly since there may be interplay between the area
conservation of the sub-cortex and the imposed steric effects at a radial
boundary. More specifically, if the lobes are not allowed to grow radially,
they may grow tangentially, making way for sharper folds in the cerebellum. To
check for this, we model the steric effects on the developing cerebellum by
incorporating a logistic function (1+tanh($x$))/2 into the energy functional,

$\begin{split}E\bigg{[}r,t,\frac{dt}{d\theta}\bigg{]}=\int
d\theta\bigg{\\{}&k_{r}(r-r_{0})^{2}-k_{t}(t-t_{0})^{2}\\\
&+\beta\bigg{(}\frac{dt}{d\theta}\bigg{)}^{2}+\frac{K_{c}}{2}(\text{tanh}[q(r-r_{c})])\bigg{\\}}.\end{split}$
(23)

The constants of the logistic function are omitted as they vanish in the
resulting Euler-Lagrange equations. Here, $K_{c}$ is the coupling constant,
$q^{-1}$ is the width of the step size of the tanh function and $r_{c}$ is the
radial position of the step. The strength of the steric interaction is taken
to be $K_{c}\;q^{2}$, a quantity whose dimensions matches those of $k_{r}$.
Renormalizing the parameters and variables by $q^{-1}$ avoids singularities in
the Euler-Lagrange equations in the limit of $q\rightarrow 0$. Dividing Eq. 23
by $k_{r}q^{-2}$, we have,

$\begin{split}\tilde{E}\bigg{[}\tilde{r},\tilde{t},\frac{d\tilde{t}}{d\theta}\bigg{]}=\int
d\theta\bigg{\\{}&(\tilde{r}-\tilde{r}_{0})^{2}-\frac{1}{c}(\tilde{t}-\tilde{t}_{0})^{2}\\\
&+\frac{1}{\rho\;c}\bigg{(}\frac{d\tilde{t}}{d\theta}\bigg{)}^{2}+\tilde{K}_{c}\;\text{tanh}[\tilde{r}-\tilde{r}_{c}]^{2}\bigg{\\}},\end{split}$
(24)

where $\tilde{E}:=E/(k_{r}q^{-2})$, $\tilde{r}:=qr$, $\tilde{t}:=qt$,
$\tilde{r}_{c}:=qr_{c}$ and $\tilde{K}_{c}=K_{c}q^{2}/k_{r}$. All parameters
and variables are now rendered dimensionless. Given that there is no explicit
dependence on $q$, we can be assured that it won’t show up in the Euler-
Lagrange equation either. The corresponding area conserving constraint is,

$\frac{1}{2}\int
d\theta(\tilde{r}-\tilde{t})^{2}=A_{0}q^{2}=\text{dimensionless constant}.$
(25)

The Euler-Lagrange equation is,

$\begin{split}&\\{(1-\epsilon)-\tilde{K}_{c}\;\text{sech}^{2}(\tilde{r}-\tilde{r}_{c})\;\text{tanh}(\tilde{r}-\tilde{r}_{c})\\}\left(\frac{d^{2}\tilde{r}}{d\theta^{2}}\right)\\\
&+\tilde{K}_{c}\;\text{sech}^{2}(\tilde{r}-\tilde{r}_{c})\;\\{2-3\;\text{sech}^{2}(\tilde{r}-\tilde{r}_{c})\\}\left(\frac{d\tilde{r}}{d\theta}\right)^{2}\\\
&+\left\\{\rho\epsilon^{2}c+(1-\epsilon)\gamma^{2}\right\\}\;\tilde{r}\\\
&+\tilde{K}_{c}\gamma^{2}\;\text{sech}^{2}(\tilde{r}-\tilde{r}_{c})\\\
&+\rho\epsilon\tilde{t_{0}}-\gamma^{2}\tilde{r}_{0}=0.\end{split}$ (26)

Here $\gamma=\sqrt{\rho(\epsilon c+1)}$. We, again, observe the
$(dr/d\theta)^{2}$ term. However, its coefficient is zero for all $r\neq
r_{c}$, thus effectively localizing its effect. Results of numerical
integration in Fig. 4a,b shows local flattening of the crests/gyri at contact
with the confining wall. Interestingly, there are no nonlocal effects of the
confining wall on the lobes, i.e., steric effects do not contribute to the
sharpness of the troughs/sulci, at least for the parameters we study. The area
conservation on the sub-cortex is not a strong enough constraint to effect the
shape of the lobes near the troughs as the sub-cortex can still change its
shape.

##  V A branching hierarchy

Figure 5: Branching morphogenesis in the cerebellum and length-scale
invariance in BWBM foliation. a) Midline saggital cross-section of mouse
cerebellum 3 days post birth. Scale bar is 200 $\mu$m. Reprinted from Fig. 7b
of Ref. Lawton  _et al._ (2019) published under the terms of the Creative
Commons Attribution. b) Idealized geometrical example of hierarchy manifests
as a fractal. c) Representation of branching morphogenesis implemented using
the geometrical idea of (b) and the 4 lobe BWBM systems with elliptical
$r_{0}$ of (d). Each lobe becomes a sub-system of its own and spawns its own
folds. Here, the first and each of the second generation systems are generated
numerically with the second generation overlaying the first after a rescaling
factor of 0.28. d,e) Folds in (d) develop transversely to the major axis of
the elliptical $r_{0}$ (blue) in contrast to (e) which has a circular $r_{0}$.
Parameters:
$c=0.0026,\epsilon=0.6,\rho=15.6,\tilde{t}_{0}=0.58,\tilde{k}_{r_{1}}=\tilde{k}_{r_{2}}$
and eccentricity $e=0,06$ respectively for the purple and green 4 lobed
systems respectively. The initial condition at $t=0$ is
[$\tilde{r},\frac{d\tilde{r}}{d\theta}]=[0.76,0]$. The parameters used for the
6 lobed system is the same as the linear BWBM values in Fig. 2. The 4 lobed
second generation lobes in (c) are generated using the parameter values
corresponding to $e=0.6$ in (b).

Mammalian cerebellums are seen to develop folds irrespective of the size of
the organ Larsell (1967). This is in contrast to small mammalian cerebrums
which do not develop folds Manyuhina  _et al._ (2014); Tallinen  _et al._
(2014). A feature of the linear BWBM model is that it offers an explanation
for the length scale independence of the folding morphogenesis in the
cerebellum. In the limit of small $\epsilon$ of the linear BWBM model, the
number of primary folds $N$ scales as $\sqrt{\rho}$. Thus, in this limit, the
number of folds developed in the cerebellum is independent of the size of the
cerebellum and is determined solely by the material elastic constants. If the
material elastic constants do not dramatically differ across mammalian
species, the linear BWBM model can potentially explain the conservation in the
number of primary folds. In this section, we discuss how the linear BWBM can
be used to describe the branching hierarchy within a given cerebellum.

As the gyri/crests of the cerebellar cortex continue to grow, the
sulci/troughs sharpen and become anchored Sudarov and Joyner (2007). The
anchoring is due to a combination of the radial glial cells and the pial
basement membrane, a thin sheet of extracellular matrix (ECM) made up of
collagen, laminin, and other ECM components Halfter  _et al._ (2002). It is
known that basement membranes stick together Baumann (2014), so as the troughs
sharpen, the pial membrane from each side of the trough begins to come into
contact and stick to reinforce the anchoring. The resulting anchoring centers
delineate the petal-like projections called lobules Marzban (2017).

We hypothesize that when the anchoring centers serve as effective boundaries
between the lobules, each lobe becomes its own subsystem. This subsystem then
consists of its own subcortex/subcore, all within the encompassing primary
cortex/core geometry. Some of these featureless subsystems go on to develop
folds of their own to, in turn, generate another generation of subsystems and
so on. See Fig. 5a. In other words, as the cerebellum continues to develop, a
branching hierarchy of subsystems emerges. This branching hierarchy is yet
another distinguishing feature of the cerebellum that sets it apart from the
cerebrum. It is to be noted that even though cerebellums of all sizes
demonstrate folds, smaller cerebellums have fewer hierarchial branchings.

The hierarchical generation of lobules within lobules points to a scale-
invariant branching process. Given that within the framework of BWBM, the
formation of crests and troughs do not depend on system size, this framework
offers a natural description for the hierarchy. As the size of the subsystems
decreases with each successive generation, crests and troughs can still form
as long as the material properties do not change. The validity of the two
dimensional model to describe growth in three dimensions remains valid since
the cerebellum retains its cylindrical symmetry during development Lawton  _et
al._ (2019).

As an idealized, purely geometric example of the hierarchy, we consider an
initial zeroth generation circle. Along its perimeter, six first generation
circles are generated. This process can proceed ad infinitum, to generate a
fractal structure with a fractal dimension of $\text{log}(3)/\text{log}(2)$.
We show four generations of this hierarchy in Fig. 5b.

In the same vein, using the linear BWBM model, we represent the hierarchy of
folds in the cerebellum in Fig. 5c. The ‘zeroth’ generation is the circular
$r_{0}$ (not shown in figure) which generates six first generation foliations.
Each lobe formed within consecutive folds is now considered an independent
subsystem and sets the length scale for the $r_{0}$ of the following, second
generation of lobes.

Fig. 5a suggests the free part of the first generation lobe i.e. the part that
is not sticking to its neighbors, to be ‘elliptical’ with the major axis of
the ellipse being parallel to the outermost exposed edge of the lobe. This is
especially evident for the L678 lobe. To accomodate this visual observation in
generating second generation folds, we employ a geometric non-linearity in the
form of an elliptical $r_{0}$ with an eccentricity $e$ within the BWBM model.
We assume there are no residual stresses and generate a 4 lobed BWBM system as
in Fig. 5d. We then use only the top half of this solution to represent the
second generation lobes in Fig. 5c.

It is also interesting to note that in Fig. 5d, the most prominent fold occurs
transversely to the horizontal major axis of the elliptical $r_{0}$. This is
simply the consequence of the system trying to minimize its energy
contribution from the radial glial springs. For comparison, we show in Fig.
5e, a four lobed system generated formed from a circular $r_{0}$.

##  VI Discussion

Inspired by cerebellar shape development, we study the effects of nonlinear
elasticity, steric confinement, and a branching hierarchy within the BWBM
model. Exploring the effects of nonlinear elasticity of the fibrous radial
glial cells, the interplay between geometry and nonlinearity is seen to give
rise to troughs sharper than the troughs obtained in the linear BWBM model and
we arrive at an asymmetry between the crests (gyri) and the troughs (sulci).
This asymmetry can be understood thus: the relatively slow growth of the sub-
cortex is taken into account by demanding area conservation of the sub-cortex
in the two-dimensional BWBM model. The associated Lagrange multiplier couples
the radius of the cerebellum and the thickness of the cortex. Nonlinear radial
springs, in association with this coupling, results in the robust quadratic
nonlinear term of the form $(dr/d\theta)^{2}$ in the shape equation. To
illustrate the role of the quadratic nonlinear term in sharpening the sulci,
we study the simple harmonic oscillator with the same form of nonlinearity and
observe its sufficiency in achieving sharp sulci. Several other nonlinearities
emerge in the shape equation including a spatially-varying effective ‘mass’
coefficient.

The perspective of cerebellum foliation as the action of a nonlinear
oscillator can be a useful one given the extensive theoretical studies of such
oscillators Nayfeh and Mook (2008); Hayashi (2014). For BWBM of the
cerebellum, the linear model with constant $r_{0}$ maps to a forced harmonic
oscillator and, for small eccentricities of $r_{0}$, maps to an unconventional
Duffing oscillator. For nonlinear $\tilde{K}_{r}(r)$, we attempt to understand
the corresponding nonlinearity in the context of the assisting-dampening
oscillator. We hope the study of cerebellar foliation as a nonlinear
oscillator problem continues to be fruitful. In a related work, the existence
of a new morphological instability in confined nonlinear elastic sheets was
found in the context of a period-doubling bifurcation, exhibiting an analogy
with parametric resonance in another nonlinear oscillator Brau  _et al._
(2011).

The period-doubling hierarchy found in the Ref. Brau  _et al._ (2011) is very
different from the new hierarchy found here. The hierarchy found here is one
due to boundary conditions in the form of anchoring centers to create sub-
regions, or sub-systems, from which the same type of scale-invariant foliation
emerges, at least in the limit of small $\epsilon$. Had the foliation
mechanism not been scale-invariant in any limit, such as with a purely elastic
system, the smaller sub-lobes would soon become featureless as the number of
foliations depend linearly on the perimeter of the sub-system. Within BWBM,
therefore, we have identified a new scale-invariant branching morphogenesis
mechanism. It is not yet clear how generic this new branching hierarchy
mechanism is in terms of moving beyond the cerebellum. Ref. Engstrom  _et al._
(2018) addressed potential applications of BWBM to two-dimensional brain
organoids Karzbrun  _et al._ (2018) and the developing retina Kolb (2018).

Referring again to Fig. 5(a), one of the second generation sublobes labelled
L45 does not branch. Perhaps the material properties are altered in this
sublobe so that features do not form. For sublobe L678, the two new sub-
sublobes are not similar in size. This could be due to changes in curvature of
the confinement from growing, surrounding tissue. So far, we have only
addressed steric, static confinement. Certainly, such variabilities from
sublobe to sublobe can be explored in less idealized conditions.

We note that within the context of purely elasticity theory, an explanation
for the thicker-sulci/thinner-gyri of the developing cerebellar cortex was
recently achieved by the addition of surface tension Riccobelli and Bevilacqua
(2020). While more energetic contributions can certainly be added to either
the purely elastic model or to the BWBM model, the recent experimental
observations needs to be incorporated in the modeling – the cells in the
cortex are motile with cellular rearrangements on the time scale of minutes
Lawton  _et al._ (2019) and the cerebellum is under tension as it develops (as
opposed to compression). These observations render a purely elastic model
suspect. However, the differential growth between the cortex and sub-cortex
remains central in both classes of models. To make progress, we need further
experimental falsification tests to rule out classes of models.

Finally, given our focus on the cerebellum here, one is led to wonder whether
some form of BWBM is applicable to the cerebrum. As mentioned previously, the
cerebrum and the cerebellum have rather different morphologies. In terms of
development in the cerebellum, the predominant growth of the cells is in the
cortex Sudarov and Joyner (2007). Many such cells migrate inward to become
part of the core of the cerebellum. In the cerebrum, much of the cell
proliferation is in the core and the progenitor cells migrate outward to
become part of the cortex Rakic (2009); Molnár  _et al._ (2019) . In this
sense, the two organs are inverse to each other. Given the presence of motile
cells in the developing cerebrum, one may wonder whether a purely elastic
approach to the developing brain is reasonable. Without a doubt, the time
scales of cell migration decide the contest between elasticity and fluidity.
Under the framework of liquid crystals, earlier work on the developing
cerebrum arrived at a prediction for the bending modulus of the cortex
Manyuhina  _et al._ (2014). This approach was based on a revised view of the
axonal tension model for the developing cerebrum Van Essen (1997, 2020). These
novel approaches have the potential to build new inroads in quantitative
biology.

MCG acknowledges useful discussions with Manu Mannattil. JMS acknowledges
financial support from NSF-DMR-1832002 and an Isaac Newton Award from the DoD.

## References

  * Ben Amar  _et al._ (2019) M. Ben Amar, P. Nassoy, and L. LeGoff, Philosophical Transactions of the Royal Society A 377, 20180070 (2019).
  * Nelson (2016) C. M. Nelson, Journal of Biomechanical Engineering 138 (2016).
  * Budday  _et al._ (2015) S. Budday, P. Steinmann III, and E. Kuhl, Frontiers in Cellular Neuroscience 9, 257 (2015).
  * Richman  _et al._ (1975) D. P. Richman, R. M. Stewart, J. W. Hutchinson, and V. S. Caviness, Science 189, 18 (1975).
  * Raghavan  _et al._ (1997) R. Raghavan, W. Lawton, S. Ranjan, and R. Viswanathan, Journal of Theoretical Biology 187, 285 (1997).
  * Bayly  _et al._ (2013) P. Bayly, R. Okamoto, G. Xu, Y. Shi, and L. Taber, Physical Biology 10, 016005 (2013).
  * Tallinen  _et al._ (2014) T. Tallinen, J. Y. Chung, J. S. Biggins, and L. Mahadevan, Proceedings of the National Academy of Sciences 111, 12667 (2014).
  * Tallinen  _et al._ (2016) T. Tallinen, J. Y. Chung, F. Rousseau, N. Girard, J. Lefèvre, and L. Mahadevan, Nature Physics 12, 588 (2016).
  * Hannezo  _et al._ (2011) E. Hannezo, J. Prost, and J.-F. Joanny, Physical Review Letters 107, 078104 (2011).
  * Shyer  _et al._ (2013) A. E. Shyer, T. Tallinen, N. L. Nerurkar, Z. Wei, E. S. Gil, D. L. Kaplan, C. J. Tabin, and L. Mahadevan, Science 342, 212 (2013).
  * Wiggs  _et al._ (1997) B. R. Wiggs, C. A. Hrousis, J. M. Drazen, and R. D. Kamm, Journal of Applied Physiology 83, 1814 (1997).
  * Li  _et al._ (2011) B. Li, Y.-P. Cao, X.-Q. Feng, and H. Gao, Journal of the Mechanics and Physics of Solids 59, 758 (2011).
  * Osborn (2008) J. W. Osborn, Journal of Theoretical Biology 255, 338 (2008).
  * Shyer  _et al._ (2017) A. E. Shyer, A. R. Rodrigues, G. G. Schroeder, E. Kassianidou, S. Kumar, and R. M. Harland, Science 357, 811 (2017).
  * Lejeune  _et al._ (2016) E. Lejeune, A. Javili, J. Weickenmeier, E. Kuhl, and C. Linder, Soft Matter 12, 5613 (2016).
  * Lejeune  _et al._ (2019) E. Lejeune, B. Dortdivanlioglu, E. Kuhl, and C. Linder, Soft Matter 15, 2204 (2019).
  * Mongera and et al. (2018) A. Mongera and et al., Nature 561, 401 (2018).
  * Jain and et al (2020) A. Jain and et al, Nature Communications 11, 5604 (2020).
  * Engstrom  _et al._ (2018) T. Engstrom, T. Zhang, A. Lawton, A. Joyner, and J. M. Schwarz, Physical Review X 8, 041053 (2018).
  * Larsell (1967) O. Larsell,  _The Comparative Anatomy and Histology of the Cerebellum_ (University of Minnesota Press, 1967).
  * Herculano-Houzel (2009) S. Herculano-Houzel, Front. Hum. Neurosci. 3, 31 (2009).
  * Lawton  _et al._ (2019) A. K. Lawton, T. Engstrom, D. Rohrbach, M. Omura, D. H. Turnbull, J. Mamou, T. Zhang, J. M. Schwarz, and A. L. Joyner, Elife 8, e45019 (2019).
  * Fernandez  _et al._ (2006) P. Fernandez, P. A. Pullarkar, and A. Ott, Biophys. J. 90, 3796 (2006).
  * Varner and Nelson (2014) V. D. Varner and C. M. Nelson, Development 141, 2750 (2014).
  * Hannezo  _et al._ (2017) E. Hannezo, C. L. G. J. Scheele, M. Moad, N. Drogo, R. Heer, R. V. Sampongna, J. van Rheenen, and B. D. Simons, Cell 171, 242 (2017).
  * Lagrange  _et al._ (2016) R. Lagrange, F. L. Jiménez, D. Terwagne, M. Brojan, and P. Reis, Journal of the Mechanics and Physics of Solids 89, 77 (2016).
  * Marzban (2017) H. Marzban,  _Development of the Cerebellum from Molecular Aspects to Diseases_ (Springer, 2017).
  * Storm  _et al._ (2005) C. Storm, J. J. Pastore, F. C. MacKintosh, T. C. Lubensky, and J. P. A, Nature 435, 195 (2005).
  * Davis (1961) H. T. Davis,  _Introduction to Nonlinear Differential and Integral Equations_ (US Government Printing Office, 1961).
  * Ermentrout (2002) B. Ermentrout,  _Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students_ (SIAM, 2002).
  * Goldstein  _et al._ (2002) H. Goldstein, C. Poole, and J. Safko,  _Classical mechanics_ (Addison Wesley, 2002).
  * Strogatz (2018) S. H. Strogatz,  _Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering_ (CRC Press, 2018).
  * Jordan  _et al._ (2007) D. Jordan, P. Smith, P. Smith, _et al._ ,  _Nonlinear Ordinary Differential Equations: An Introduction for Scientists and engineers_ , Vol. 10 (Oxford University Press on Demand, 2007).
  * Stokes (1880) G. G. Stokes, Transactions of the Cambridge Philosophical Society (1880).
  * Nie  _et al._ (2010) J. Nie, L. Guo, G. Li, C. Faraco, L. S. Miller, and T. Liu, Journal of Theoretical Biology 264, 467 (2010).
  * Manyuhina  _et al._ (2014) O. Manyuhina, D. Mayett, and J. Schwarz, New Journal of Physics 16, 123058 (2014).
  * Sudarov and Joyner (2007) A. Sudarov and A. L. Joyner, Neural Development 2, 26 (2007).
  * Halfter  _et al._ (2002) W. Halfter, S. Dong, Y.-P. Yip, M. Willem, and U. Mayer, J. Neuroscience 22, 6029 (2002).
  * Baumann (2014) K. Baumann, Nature Reviews Molecular Cell Biology 15, 767 (2014).
  * Nayfeh and Mook (2008) A. H. Nayfeh and D. T. Mook,  _Nonlinear Oscillations_ (John Wiley & Sons, 2008).
  * Hayashi (2014) C. Hayashi,  _Nonlinear Oscillations in Physical Systems_ (Princeton University Press, 2014).
  * Brau  _et al._ (2011) F. Brau, H. Vandeparre, A. Sabbah, C. Poulard, A. Boudaoud, and P. Damman, Nature Physics 7, 56 (2011).
  * Karzbrun  _et al._ (2018) E. Karzbrun, A. Kshirsagar, S. R. Cohen, J. H. Hanna, and O. Reiner, Nature Physics 14, 515 (2018).
  * Kolb (2018) H. Kolb, in  _Webvision: The Organization of the Retina and Visual System_ , edited by H. Kolb, R. Nelson, E. Fernandez, and B. Jones (Webvision, 2018) Chap. 1.
  * Riccobelli and Bevilacqua (2020) D. Riccobelli and G. Bevilacqua, Journal of the Mechanics and Physics of Solids 134, 103745 (2020).
  * Rakic (2009) P. Rakic, Nature Reviews Neuroscience 10, 724 (2009).
  * Molnár  _et al._ (2019) Z. Molnár, G. J. Clowry, N. Šestan, A. Alzu’bi, T. Bakken, R. F. Hevner, P. S. Hüppi, I. Kostović, P. Rakic, E. Anton, _et al._ , Journal of Anatomy 235, 432 (2019).
  * Van Essen (1997) D. C. Van Essen, Nature 385, 313 (1997).
  * Van Essen (2020) D. C. Van Essen, Proceedings of the National Academy of Sciences 117, 32868 (2020).

##  VII Appendix

###  VII.1 Uncoupling general Euler-Lagrange
equations for $r,t$

The uncoupled nonlinear differential equation in $r$ for the case of having
$k_{r}:=k_{r}(r)$ is shown in Eq. 14. The form of such an equation is,
however, dependent on the choice of nonlinear coefficients. For general
coupled Euler-Lagrange equations, we have

$f(r,t)=0,$ (27) $t^{\prime\prime}=g(r,t).$ (28)

Here the prime superscript indicates a derivative with respect to $\theta$.
The derivative of Eq. 27 with respect to $\theta$ gives,

$\frac{\partial f(r,t)}{\partial r}r^{\prime}+\frac{\partial f(r,t)}{\partial
t}t^{\prime}=0.$ (29)

Another derivative yields,

$\begin{split}\frac{\partial f(r,t)}{\partial
r}r^{\prime\prime}+\frac{\partial f(r,t)}{\partial
t}t^{\prime\prime}+\frac{\partial^{2}f(r,t)}{\partial r^{2}}r^{\prime
2}&+2\;\frac{\partial^{2}f(r,t)}{\partial r\partial t}r^{\prime}t^{\prime}\\\
&+\frac{\partial^{2}f(r,t)}{\partial t^{2}}t^{\prime 2}=0.\end{split}$ (30)

Eq. 27,29 can be solved to find
$t:=t(r),\;t^{\prime}=t^{\prime}(r,r^{\prime})$. The uncoupled differential
equation in $r$ would then be,

$\begin{split}&\frac{\partial f(r,t(r))}{\partial
r}r^{\prime\prime}+\frac{\partial f(r,t(r))}{\partial
t}t^{\prime\prime}(r)+\frac{\partial^{2}f(r,t(r))}{\partial r^{2}}r^{\prime
2}\\\ &+2\;\frac{\partial^{2}f(r,t(r))}{\partial r\partial
t}r^{\prime}t^{\prime}(r)+\frac{\partial^{2}f(r,t(r))}{\partial
t^{2}}t^{\prime 2}(r)=0.\end{split}$ (31)

The last two terms for Eq. 14 is absent and this particular case of nonlinear
differential equation is of the form,

$\frac{\partial f(r,t(r))}{\partial
r}r^{\prime\prime}+\frac{\partial^{2}f(r,t(r))}{\partial r^{2}}r^{\prime
2}+\frac{\partial f(r,t(r))}{\partial t}g(r,t(r))=0,$ (32)

where the implicit dependence - $t(r)$ is to be taken into account after
evaluating the partial derivatives in $t,r$. This renders an uncoupled
differential equation of motion for $r$.

###  VII.2 Higher order corrections to the Simple Harmonic Oscillator

With no loss of generality, we can study the simple harmonic oscillator with
cubic and quartic corrections to the potential energy as,

Figure 6: Lindstedt-Poincaré perturbation method for the cubic and quartic
energy corrections to the SHO. $\mathcal{O}(\epsilon^{2})$ curve follows the
leading edge of the numerical solution but doesn’t capture the offset in
troughs and the relative difference in widths of crests and troughs. Here
$\epsilon=0.6,\Gamma^{\prime}=0.25$.

$\frac{d^{2}x}{dt^{2}}=-x-x^{2}-\Gamma^{\prime}x^{3},$ (33)

a single parameter equation where all variables are nondimensionalized. To
evolve a perturbative scheme, we treat the last two terms of the above
equation as a perturbation to the simple harmonic oscillator. To that effect,
with $\epsilon>0$ we have,

$\frac{d^{2}x}{dt^{2}}=-x-\epsilon\;(x^{2}-\Gamma^{\prime}x^{3}).$ (34)

The Lindstedt-Poincaré method Jordan  _et al._ (2007) introduces angular
frequency $\omega$ by change of variable $\tau=\omega t$. This allows for
coupling between the frequency and amplitude. Eq. 34 is now written as,

$\omega^{2}\frac{d^{2}x}{d\tau^{2}}+x+\epsilon
x^{2}-\epsilon\Gamma^{\prime}x^{3}=0.$ (35)

The series expansions for $x$ and $\omega$ are,

Figure 7: Phase portraits of BWBM and ADO systems. Blue lines are position
nullclines: $dr/dt=0$ and red lines are velocity nullclines:
$d^{2}r/dt^{2}=0$. The tangent to the flow lines are vertical when they pass
the position nullcline and horizontal when they pass the velocity nullcline.
a)Linear BWBM b) ADO c) Cubic and quartic energy corrections to the SHO d)
BWBM with nonharmonic springs. In (b,c,d) left-right asymmetric orbits where
there is a steeper fall on the left translates to sharper troughs in position-
time space. We see left-right asymmetric orbits with or without bent
nullclines (see b,d & c).

$\begin{split}x(t)&=x_{0}(t)+\epsilon x_{1}(t)+\epsilon^{2}x_{2}(t)+...\\\
\omega&=1+\epsilon\omega_{1}+\epsilon^{2}\omega_{2}+...\end{split}$ (36)

The boundary conditions are,

$\begin{split}x_{0}(0)&=A_{0},\dot{x}_{0}(0)=0,\\\
x_{i}(0)&=0,\dot{x}_{i}(0)=0~{}~{}i>0.\end{split}$ (37)

$\mathcal{O}(\epsilon^{0})$
We substitute the series expansions of Eq. 36 in Eq. 35. All the terms at
every given order of $\epsilon$ must add up to zero for Eq. 35 to hold. The
response equation at the zeroth order of $\epsilon$ then is,

$\frac{d^{2}x_{0}}{d\tau^{2}}+x_{0}=0,$ (38)

which is just the equation of motion of the simple harmonic oscillator. The
zeroth order response that obeys the boundary conditions 37 is,

$x_{0}(\tau)=A_{0}\text{cos}(\tau).$ (39)

$\mathcal{O}(\epsilon^{1})$
At the first order of $\epsilon$ we have,

$\frac{d^{2}x_{1}}{d\tau^{2}}+x_{1}=x_{0}^{2}+\Gamma^{\prime}x_{0}^{3}-2\omega_{1}\frac{d^{2}x_{0}}{d\tau^{2}}.$
(40)

Solutions of lower order response equations are used to find solutions of
higher order response equations. $\omega_{1}$ is fixed by the demand that the
solution to Eq. 40 be periodic. For this, upon substituting the zeroth order
solution (see Eq. 39) in Eq. 40 we set the coefficient of cos$(\tau)$ in the
resulting RHS of Eq. 40 to vanish. We obtain

$\omega_{1}=\frac{-3\Gamma^{\prime}a_{0}^{2}}{8}.$ (41)

The amplitude-frequency dependence can be seen here. The first order
correction to $x$ which obeys the boundary conditions is,

$\begin{split}x_{1}(\tau)=A_{0}^{2}\left(\frac{\Gamma^{\prime}A_{0}}{32}-\frac{1}{3}\text{cos}(\tau)\right)+\frac{A_{0}^{2}}{2}&-\frac{A_{0}^{2}}{6}\text{cos}(2\tau)\\\
&-\frac{\Gamma^{\prime}A_{0}^{3}}{32}\text{cos}(3\tau).\end{split}$ (42)

$\mathcal{O}(\epsilon^{2})$
At the second order of $\epsilon$, we have,

$\frac{d^{2}x_{2}}{d\tau^{2}}+x_{2}=3\Gamma^{\prime}x_{0}^{2}x_{1}+2x_{0}x_{1}-(\omega_{1}^{2}+2\omega_{2})\frac{d^{2}x_{0}}{d\tau^{2}}-2\omega_{1}\frac{d^{2}x_{1}}{d\tau^{2}}.$
(43)

Carrying out the same procedure at this order, we have,

$\begin{split}\omega_{2}=\frac{1}{384}\bigg{(}-9a_{0}^{4}\Gamma^{\prime
2}+144a_{0}^{3}\Gamma^{\prime}&-12a_{0}^{2}\Gamma^{\prime}\omega_{1}-160a_{0}^{2}\\\
&+128a_{0}\omega_{1}-192\omega_{1}^{2}\bigg{)},\end{split}$ (44)

and

$\begin{split}x_{2}(\tau)=&-\frac{1}{384}a_{0}^{2}[-64\text{cos}(2\tau)(3a_{0}^{2}\Gamma^{\prime}-2a_{0}-8\omega_{1})\\\
&+a_{0}(\text{cos}(3\tau)(9a_{0}^{2}\Gamma^{\prime
2}+96a_{0}\Gamma^{\prime}+216\Gamma^{\prime}\omega_{1}+64)\\\
&+60a_{0}\Gamma^{\prime}\text{cos}(4\tau)+9a_{0}^{2}\Gamma^{\prime
2}\text{cos}(5\tau)-252a_{0}\Gamma^{\prime}\\\ &+128)].\end{split}$ (45)

The perturbative expansions are compared with numerical solution in Fig. 6.

###  VII.3 Phase portraits

Closed orbits in phase-space translate to periodic motion in the position-time
space. Plotting orbits in phase space for different initial conditions builds
the phase portrait of the system. The $x\;(dx/dt)$ nullcline are the locus of
points where $dx/dt=0$ ($d^{2}x/dt^{2}=0$). The nullclines are generally used
to divide the phase-portrait into regions where the tangent of the orbit
points in the same general direction (NW, NE, SE or SW). The orbit has a
vertical (horizontal) tangent when it passes through the $x$ ($dx/d\theta$)
nullcline. The intersection of the nullclines gives the equilibrium point
where both the ‘velocity’ and ‘acceleration’ of the system is zero. The time-
reversal symmetry of the governing differential equations (see Eq. 14, 15, 26)
guarantees closed orbits at points close enough to the equilibrium point. This
is verified in the phase-portraits in Fig. 7.