# Higher-curvature generalization of Eguchi-Hanson spaces

Cristóbal Corral<EMAIL_ADDRESS>Instituto de Ciencias Exactas y Naturales,
Facultad de Ciencias, Universidad Arturo Prat, Avenida Arturo Prat Chacón
2120, 1110939, Iquique, Chile  Daniel Flores-Alfonso
<EMAIL_ADDRESS>Departamento de Física, CINVESTAV-IPN, A.P. 14-740,
C.P. 07000, Ciudad de México, Mexico  Gastón Giribet<EMAIL_ADDRESS>Department of Physics, New York University, 726 Broadway, New York, NY10003,
USA.  Physics Department, University of Buenos Aires and IFIBA-CONICET, Ciudad
Universitaria, Pabellón 1 (1428), Buenos Aires, Argentina  Julio Oliva
<EMAIL_ADDRESS>Departamento de Física, Universidad de Concepción, Casilla,
160-C, Concepción, Chile.

###### Abstract

We construct higher-dimensional generalizations of the Eguchi-Hanson
gravitational instanton in the presence of higher-curvature deformations of
general relativity. These spaces are solutions to Einstein gravity
supplemented with the dimensional extension of the quadratic Chern-Gauss-
Bonnet invariant in arbitrary even dimension $D=2m\geq 4$, and they are
constructed out of non-trivial fibrations over $(2m-2)$-dimensional Kähler-
Einstein manifolds. Different aspects of these solutions are analyzed; among
them, the regularization of the on-shell Euclidean action by means of the
addition of topological invariants. We also consider higher-curvature
corrections to the gravity action that are cubic in the Riemann tensor and
explicitly construct Eguchi-Hanson type solutions for such.

##  I Introduction

Discovered in 1979 by Tohru Eguchi and Andrew Hanson Eguchi and Hanson (1978,
1979a), and independently by Eugenio Calabi Calabi (1979) the same year, the
Eguchi-Hanson metric represents an interesting example of gravitational
instanton Hawking (1977); Eguchi and Freund (1976); Gibbons and Hawking
(1978); Page (1978); Hawking and Pope (1978); Gibbons and Hawking (1979);
Eguchi and Hanson (1979b); Eguchi  _et al._ (1980): It is a Ricci flat metric,
and it corresponds to a non-compact, self-dual, Euclidean solution to Einstein
equations in vacuum and with vanishing cosmological constant ($\Lambda=0$);
although an analog of the solution with non-vanishing $\Lambda$ also exists
Pedersen (1985). When $\Lambda=0$, the space is asymptotically locally
Euclidean (ALE). More precisely, it is asymptotically
$\mathbb{R}^{4}/\mathbb{Z}_{2}$. It admits a non-singular case whose topology
is the cotangent bundle of the 2-sphere
$\text{T}^{*}(S^{2})=\text{T}^{*}(\mathbb{CP}^{1})$, being $\mathbb{R}\times
S^{2}$ at the so-called origin ($r=a$ in the coordinates (1)-(3) below with
$m=2$ and $\Lambda=0$) and asymptotically $\mathbb{RP}^{3}$ (at $r=\infty$).
The projective space $\mathbb{RP}^{3}=S^{3}/\mathbb{Z}_{2}$ is isomorphic to
the 3-dimensional rotation group $SO(3)\simeq SU(2)/\mathbb{Z}_{2}$, and the
isometry group of the Eguchi-Hanson space is thus given by $U(1)\times
SU(2)/\mathbb{Z}_{2}\simeq U(2)$. As it is for the $K3$ Calabi-Yau space, the
reduced holonomy group of Eguchi-Hanson space is $SU(2)$, and this makes it
possible to approximate the former space by combinations of copies of the
latter. This is related to the fact that Eguchi-Hanson may be regarded as the
resolution of the $A_{1}$ singularity, according to the ADE classification. In
fact, Eguchi-Hanson has interesting applications in string theory in relation
to the resolutions of orbifold singularities. It also has other interesting
applications, such as serving as a seed for constructing higher-dimensional
spaces with special holonomy Higashijima  _et al._ (2002a, b, c); Clarkson and
Mann (2006a, b); Giribet and Santillan (2007); Wong and Mann (2012).
Multicenter, asymptotically $\mathbb{R}^{4}/\mathbb{Z}_{k}$ generalization of
the Eguchi-Hanson instanton can also be constructed Gibbons and Hawking
(1978).

In this paper, we are interested in studying the higher-dimensional extension
of the Eguchi-Hanson space, given by the metric (1)-(3) below. Considering the
higher-dimensional generalization of Einstein spaces, however, makes it
natural to consider, in addition to the Einstein tensor, higher-curvature
contributions yielding second-order, covariantly conserved, rank-2 tensors in
the field equations, cf. Lovelock (1971). Such contributions do not exist in 4
dimensions Lanczos (1938) but are certainly possible in dimension 5 and higher
Lovelock (1972). Here, we will take them into account. Higher-dimensional,
higher-curvature gravitational instantons have been studied before; see for
instance Hendi and Dehghani (2008); Hendi  _et al._ (2012); Bueno  _et al._
(2018); Corral  _et al._ (2019, 2021) and references therein and thereof.
Here, we will explicitly construct generalizations of Eguchi-Hanson space in
arbitrary even dimension $D$ and with quadratic curvature contributions. We
will also consider examples with cubic curvature corrections.

The paper will be organized as follows: In Sec. II, we will study the higher-
dimensional analog of the Eguchi-Hanson space with cosmological constant. In
Sec. III, we will introduce the higher-curvature terms we will consider. In
Sec. IV, we will present a six-dimensional version of the Eguchi-Hanson space,
including the higher-curvature modification. We will discuss different version
of the solution, exhibiting diverse topologies in the base manifold. In Sec.
V, we will discuss the topological invariants and the regularized on-shell
action associated to the six-dimensional solution. The $D$-dimensional
solution will be presented in Sec. VI, where, in addition to the quadratic
corrections; we will also consider cubic terms in the Riemann tensor in Sec.
VII.

##  II Higher-dimensional Eguchi-Hanson in Einstein gravity

For the sake of completeness, let us review first the $D$-dimensional Eguchi-
Hanson solution in Einstein gravity. Inhomogeneous Einstein metrics were
constructed as the non-trivial fibration over arbitrary Kähler-Einstein
manifolds in Ref. Hawking and Pope (1978) (see also Page and Pope (1987)) in
even higher dimensions $D=2m$ with $m\in\mathbb{Z}$. A particular case is the
higher-dimensional generalization of the Eguchi-Hanson metric Eguchi and
Hanson (1979a), whose line element is given by

$\displaystyle\text{d}s^{2}=\frac{r^{2}}{4}f(r)\left(\text{d}\tau+\mathcal{B}\right)^{2}+\frac{\text{d}r^{2}}{f(r)}+\frac{r^{2}}{4}\text{d}\Sigma^{2}\,,$
(1)

where $\mathcal{B}=\mathcal{B}_{\mu}\text{d}x^{\mu}$ denotes the Kähler
potential 1-form such that $\Omega=\text{d}\mathcal{B}$ defines the symplectic
real form associated to the ($2m-2$)-dimensional Kähler base manifold with
line element $\text{d}\Sigma^{2}$. In four dimensions, when the base manifold
is $\mathbb{S}^{2}$, the metric (1) can be written in terms of the left-
invariant Maurer-Cartan forms of $SU(2)$. In that case, the isometry group is
$U(1)\times SU(2)/\mathbb{Z}_{2}$ and this solution is usually referred to as
the gravitational analog of the Belavin-Polyakov-Schwartz-Tyupkin instanton in
Yang-Mills theory Belavin  _et al._ (1975).

For concreteness, let us consider that the base space is the complex
projective space $\mathbb{CP}^{k}$, with $k=m-1$. Then, the metric (1) is a
solution of the $2m$-dimensional Einstein field equations with cosmological
constant,

$\displaystyle R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R+\Lambda g_{\mu\nu}=0\,,$ (2)

provided that the metric function $f(r)$ is given by

$\displaystyle
f(r)=\frac{2}{m}\left[1-\left(\frac{a}{r}\right)^{2m}\right]-\frac{1}{2\left(m^{2}-1\right)}\,\Lambda
r^{2}\,,$ (3)

where $a$ is an integration constant. This metric is completely regular
provided the parameters are chosen in such a way that there is no conical
singularity at the bolt; the latter being defined as the ($2m-2$)-dimensional
set of fixed points located at the locus $f(r_{b})=0$. Since the metric is
positive definite in Euclidean signature, the range of the radial coordinate
is $r\in\mathbb{R}_{\geq|r_{b}|}$. The regularity condition can be achieved by
demanding periodicity on the $\tau$ coordinate, namely
$\tau\sim\tau+\beta_{\tau}$, with $\beta_{\tau}$ defined through

$\displaystyle\beta_{\tau}=\frac{8\pi}{r_{b}\,f^{\prime}|_{r=r_{b}}}\,,$ (4)

where the prime denotes differentiation with respect to the radial coordinate
$r$.

In four dimensions ($m=2$) with a base manifold of $\mathbb{S}^{2}$ topology,
the metric (1) is endowed with a Misner string Misner (1963) along the
$z$-axis as noticed in Ref. Hawking and Hunter (1999). However, its position
can be made unobservable by demanding extra conditions on the periodicity of
$\tau$ Misner (1963). For instance, the topology at $r=$ constant
hypersurfaces becomes $\mathbb{S}^{3}/\mathbb{Z}_{n}$ if $0\leq\tau\leq
4\pi/n$. Thus, the absence of conical singularities alongside the
unobservability of the Misner string demands Pedersen (1985)

$\displaystyle\frac{8\pi}{r_{b}\,f^{\prime}|_{r=r_{b}}}=\frac{4\pi}{n}\;\;\;\to\;\;\;r_{b}=\sqrt{-\frac{2(n-2)}{\Lambda}}\,.$
(5)

The case with $\Lambda=0$ leads to the original Eguchi-Hanson metric where
$\beta_{\tau}=2\pi$ and $r\geq a$, representing a space with
$\mathbb{S}^{3}/\mathbb{Z}_{2}$ topology. However, depending on whether the
cosmological constant is positive or negative, different topologies can be
obtained. As shown in in Ref. Chen and Zhang (2021), the value of $r_{b}$ in
Eq. (5) is modified in non-Einstein spaces with negative constant scalar
curvature. Indeed, this is also the case in presence of higher-curvature
corrections as we study here.

In four dimensions, this metric is globally (anti-)self dual in absence of the
cosmological constant ($\Lambda=0$). In that case, $\beta_{\tau}=2\pi$ and the
action and topological invariants were obtained in Refs. Eguchi and Hanson
(1979a, 1978). Specifically, its Eucliden on-shell action vanishes, its Euler
characteristic is $\chi=2$, and its Hirzebruch signature is $\tau=-1$. This
implies that there exists only one anti-self-dual $2$-form with the Atiyah-
Patodi-Singer boundary conditions Atiyah  _et al._ (1975, 1976a, 1976b).
Moreover, since the index of the Dirac operator vanishes, there is no axial
asymmetry whatsoever between left and right-handed Dirac spinors.
Nevertheless, there is an excess of two negative-chiralty spin $3/2$ spinors
as it can be seen from the index of Rarita-Schinger operator Eguchi and Hanson
(1979a, 1978).

The Euclidean Taub-NUT metrics Hawking (1977); Page (1978) and the Eguchi-
Hanson metric have similar off-diagonal components; they signal non-trivial
circle bundles. Moreover, the Taub-bolt and Eguchi-Hanson spaces have the same
type of degenerate submanifolds (bolts), although their global structures are
different. Additionally, the self-dual Taub-NUT metric is hyper-Kähler, while,
that of Eguchi-Hanson is Kähler. In the presence of a negative cosmological
constant, these two spaces differ in their asymptotic behavior: while the
Taub-NUT metric is asymptotically locally hyperbolic, the Eguchi-Hanson metric
is not Chen and Zhang (2021). These features represent the main differences
between Taub-NUT and Eguchi-Hanson spaces.

In the following, we study how the higher-dimensional Eguchi-Hanson instanton
is modified by the presence of higher-curvature corrections. In particular, we
focus on terms of the Lovelock’s series such that the field equations remain
of second order.

##  III Einstein-Gauss-Bonnet gravity

On a $D$-dimensional manifold $\mathcal{M}$, the most general metric theory
that is invariant under diffeomorphisms and local Lorentz transformations,
which leads to second order field equations is known as the Lanczos-Lovelock
theory of gravity Lanczos (1938); Lovelock (1971). Its dynamics is dictated by
the action principle

$\displaystyle
I[g_{\mu\nu}]=\sum_{p=0}^{\left[\frac{D-1}{2}\right]}\int_{\mathcal{M}}\text{d}^{D}x\sqrt{|g|}\;\alpha_{p}\,\mathscr{L}^{(p)}\,,$
(6)

where $g=\det g_{\mu\nu}$ is the metric determinant, $[\dots]$ represent the
integer part, $\alpha_{p}$ are the Lovelock’s couplings, and the $p$-th order
Lagrangian is defined as

$\displaystyle\mathscr{L}^{(p)}=\frac{1}{2^{p}}\,\delta^{\mu_{1}\ldots\mu_{2p}}_{\nu_{1}\ldots\nu_{2p}}R^{\nu_{1}\nu_{2}}_{\
\ \ \mu_{1}\mu_{2}}\dots R^{\nu_{2p-1}\nu_{2p}}_{\ \ \ \mu_{2p-1}\mu_{2p}}\,,$
(7)

where
$\delta^{\mu_{1}\ldots\mu_{k}}_{\nu_{1}\ldots\nu_{k}}=k!\,\delta_{[\nu_{1}}^{[\mu_{1}}\dots\delta^{\mu_{k}]}_{\nu_{k}]}$
denotes the completely antisymmetrized Kronecker delta of rank $k$. The field
equations are obtained by performing stationary variations of the action (6)
with respect of the metric, giving

$\displaystyle
E^{\mu}_{\nu}\equiv-\sum_{p=0}^{\left[\frac{D-1}{2}\right]}\frac{\alpha_{p}}{2^{p+1}}\,\delta^{\mu\mu_{1}\ldots\mu_{2p}}_{\nu\nu_{2}\ldots\nu_{2p}}R^{\nu_{1}\nu_{2}}_{\
\ \ \mu_{1}\mu_{2}}\dots R^{\nu_{2p-1}\nu_{2p}}_{\ \ \
\mu_{2p-1}\mu_{2p}}=0\,.$ (8)

These equations of motion are at most of second order in derivatives of the
metric and they propagate $D(D-3)/2$ degrees of freedom Henneaux  _et al._
(1990) around flat space; the same as general relativity in a $D$-dimensional
space, the latter being the particular case when $\alpha_{p\geq 2}=0$. Since
we are interested in non-trivial fibrations of Kähler-Einstein manifolds, we
shall focus on even-dimensional manifolds henceforth111Generalizations of
Eguchi-Hanson spaces in odd dimensions were considered in Wong and Mann
(2012).

The so-called Einstein-Gauss-Bonnet theory can be obtained by truncating the
Lovelock series to second order in $p$. Nevertheless, we will also keep the
highest-order term in the Lovelock series Lovelock (1972) which corresponds to
the $2m$-dimensional Euler density Chern (1945); namely,
$\mathscr{L}^{(m)}=\mathcal{E}_{2m}$ (see Eq. (21) below). Being a topological
invariant, the addition of the latter to the action does not modify the bulk
dynamics, but it plays a crucial role in renormalizing the Euclidean on-shell
action as well as in the computation of conserved charges in asymptotically
locally AdS (AlAdS) spaces, providing a bulk counterterm once its
corresponding coefficient is fixed Aros  _et al._ (2000a, b); Mora  _et al._
(2004); Kofinas and Olea (2006); Olea (2007); Arenas-Henriquez  _et al._
(2019). Moreover, in four dimensions, when evaluated in configurations with
(anti-)self-dual Weyl tensor, the topologically renormalized action yields an
Euclidean on-shell action that is proportional to the Pontryagin index
Ciambelli  _et al._ (2021); Corral  _et al._ (2021), similarly as it happens
for instantons in Yang-Mills theory. For non-conformally flat boundaries,
additional counterterms based on the conformal completion of Einstein-AdS
gravity may also be needed in order to render the Euclidean on-shell action
and conserved charges finite Anastasiou  _et al._ (2020, 2021a, 2021b). Adding
the Pontryagin density into the four-dimensional action with fixed coupling is
also an interesting possibility, since it sets any the (anti-)self dual
configurations as the ground state of the theory Miskovic and Olea (2009);
Araneda  _et al._ (2016).

Setting $\alpha_{0}=-2\kappa\Lambda$, $\alpha_{1}=\kappa$,
$\alpha_{2}=\kappa\alpha$, and $\kappa=(16\pi G_{N})^{-1}$ with $G_{N}$ being
the Newton’s constant, the action becomes

$\displaystyle I_{\rm
EGB}[g]=\kappa\int_{\mathcal{M}}\text{d}^{2m}x\sqrt{|g|}\,\left(R-2\Lambda+\alpha\,\mathcal{G}+\zeta\,\mathcal{E}_{2m}\right)\,,$
(9)

where $\alpha$ and $\zeta$ are dimensionful coupling constants, and
$\mathcal{G}$ is the Chern-Gauss-Bonnet term, namely,

$\displaystyle\mathcal{G}=\frac{1}{4}\delta^{\mu\nu\lambda\rho}_{\alpha\beta\gamma\delta}R^{\alpha\beta}_{\
\mu\nu}R^{\gamma\delta}_{\
\lambda\rho}=R^{2}-4R^{\mu}_{\nu}R^{\nu}_{\mu}+R^{\mu\nu}_{\ \lambda\rho}R_{\
\mu\nu}^{\lambda\rho}\,.$ (10)

In four dimensions, this term represents the bulk piece of the Euler
characteristic [see Eq. (21) below] and, therefore, it does not contribute to
the bulk dynamics Lanczos (1938). In dimensions higher than four, in contrast,
this term does modify the field equations, leading to

$\displaystyle
R^{\mu}_{\nu}-\frac{1}{2}\delta^{\mu}_{\nu}R+\Lambda\delta^{\mu}_{\nu}+\alpha
H^{\mu}_{\nu}=0\,,$ (11)

where $H^{\mu}_{\nu}$ is the contribution to the field equations of the Chern-
Gauss-Bonnet term defined through

$\displaystyle
H^{\mu}_{\nu}=-\frac{1}{8}\delta^{\mu\mu_{1}\ldots\mu_{4}}_{\nu\nu_{1}\ldots\nu_{4}}R^{\nu_{1}\nu_{2}}_{\
\ \mu_{1}\mu_{2}}R^{\nu_{3}\nu_{4}}_{\ \ \mu_{1}\mu_{4}}=2R^{\mu\rho}_{\
\sigma\lambda}R_{\ \nu\rho}^{\sigma\lambda}-4R_{\rho}^{\sigma}R^{\mu\rho}_{\
\nu\sigma}+2RR^{\mu}_{\nu}-4R^{\mu}_{\lambda}R^{\lambda}_{\nu}-\frac{1}{2}\delta^{\mu}_{\nu}\mathcal{G}\,.$
(12)

Different analytic solutions to the field equations (11) are known; these
include black holes Boulware and Deser (1985); Cai (2002), wormhole geometries
Dotti  _et al._ (2007a, b), Taub-NUT type gravitational instantons Hendi and
Dehghani (2008) and many other examples. In the next section, we provide a
novel stationary solution to the field equations in Euclidean signature which
represents a higher-dimensional generalization of the Eguchi-Hanson metric.

##  IV Six-dimensional Eguchi-Hanson solution

Let us begin by considering the six-dimensional case, which is the simplest in
which the quadratic curvature terms contribute in deforming the geometry in
presence of an Kähler-Einstein base manifold. To solve the field equations
(11), we focus on a six-dimensional Eguchi-Hanson-inspired ansatz (1). In
particular, we consider different topologies for the latter, these being
fibrations over $\mathbb{CP}^{2}$, $\mathbb{CH}^{2}$, $\mathbb{T}^{4}$,
$\mathbb{S}^{2}\times\mathbb{S}^{2}$, and
$\mathbb{H}^{2}\times\mathbb{H}^{2}$. More details of these geometries can be
found in the next section. With these Kähler-Einstein base manifolds, the
field equations (11) turn out to be solved by the metric function

$\displaystyle
f_{\pm}(r)=\frac{2}{3}\left[\gamma+\frac{r^{2}}{16\alpha}\left(1\pm\sqrt{\left(1+3\Lambda\alpha\right)-\frac{128\sigma\alpha^{2}}{r^{4}}+\frac{32\alpha
a^{6}}{r^{8}}}\right)\right]\,,$ (13)

where $a$ is an integration constant, $\gamma$ is related to the constant
curvature of the base manifold, and $\sigma$ measures whether the latter is
conformally flat ($\sigma=0$) or not ($\sigma=1$).222This is similar to what
happens for solutions of the Einstein-Gauss-Bonnet field equations in the
simpler geometries of the form $M_{2}\times_{w}\Sigma_{d-2}$, where
$\times_{w}$ stands for a warped product. In that case, when $\Sigma_{d-2}$ is
not conformally flat, there is an extra constant in the lapse function –
analogous to $\sigma$ here– which leads to a slow decay of the metric at
infinity relative to the $r$-dependent behavior that depends on the
integration constant $a$, cf. Ref. Dotti and Gleiser (2005). Specifically, the
values of these parameters can be summarized as follows

| $\;\;\;\mathbb{CP}^{2}\;\;\;$ | $\;\;\;\mathbb{CH}^{2}\;\;\;$ | $\;\;\;\;\;\mathbb{T}^{4}\;\;\;\;\;$ | $\;\mathbb{S}^{2}\times\mathbb{S}^{2}\;$ | $\;\mathbb{H}^{2}\times\mathbb{H}^{2}\;$
---|---|---|---|---|---
$\;\;\;\gamma\;\;\;$ | $1$ | $-1$ | $0$ | $1$ | $-1$
$\;\;\;\sigma\;\;\;$ | $0$ | $0$ | $0$ | $1$ | $1$
(14)

The metric has a bolt at $r=r_{b}$ defined as the largest root of the
polynomial $f(r_{b})=0$ in Eq. (13). Additionally, since the Euclidean
signature implies that the metric is positive definite, the range of the
radial coordinate is restricted to be $r\in\mathbb{R}_{\geq|r_{b}|}$. The
absence of conical singularities is guaranteed as long as Eq. (4) holds for
the solution (13). Asymptotically, the metric behaves as

$\displaystyle f_{\pm}(r)$
$\displaystyle\simeq\frac{\left(1\pm\sqrt{1+3\Lambda\alpha}\right)r^{2}}{24\alpha}+\frac{2\gamma}{3}\mp\frac{8\sigma\alpha}{3\sqrt{1+3\Lambda\alpha}\;r^{2}}$
$\displaystyle\pm\frac{2\left[\left(1+3\Lambda\alpha\right)a^{6}-128\sigma^{2}\alpha^{3}\right]}{3\left(1+3\Lambda\alpha\right)^{3/2}\,r^{6}}+\mathcal{O}(r^{-10})\,.$
(15)

Thus, one can directly notice that
$\left(1\pm\sqrt{1+3\Lambda\alpha}\right)/(24\alpha)$ plays the role of
effective cosmological constant, as long as $\alpha\neq 0$. The solution is
neither asymptotically conformally flat nor asymptotically locally AdS. In
fact, the Weyl tensor squared behaves like

$\displaystyle W^{\mu\nu}_{\ \ \lambda\rho}W^{\lambda\rho}_{\ \
\mu\nu}\simeq\frac{2\left[2+3\Lambda\alpha\pm
2\sqrt{1+3\Lambda\alpha}\right]}{5\alpha^{2}}+\mathcal{O}\left(r^{-8}\right)\;\;\;\;\;\mbox{as}\;\;\;\;\;r\to\infty\,,$
(16)

and it does not vanish at infinity. Notice that the branch $f_{-}(r)$ of the
solution (13) is continuously connected to the solution of Einstein gravity
[cf. Eq. (3)] in the limit $\alpha\to 0$, since

$\displaystyle
f_{-}(r)\simeq\frac{2}{3}\left(\gamma-\frac{a^{6}}{r^{6}}\right)-\frac{\Lambda
r^{2}}{16}+\mathcal{O}(\alpha/r^{2})\,.$ (17)

Also notice that, on the curve of the parameter space defined by
$\Lambda\alpha=-1/3$, the solution acquires a rather simple form; namely,

$\displaystyle
f_{\pm}(r)=\frac{2}{3}\left[\gamma+\frac{r^{2}}{16\alpha}\left(1\pm\frac{1}{3\Lambda
r^{2}}\sqrt{-\frac{96\Lambda a^{6}}{r^{4}}-128\sigma}\right)\right]\,.$ (18)

This is the confluent point where $f_{+}(r)=f_{-}(r)$ when $a=\sigma=0$ and,
in many aspects, it exhibits certain features analogous to the so-called
Chern-Simons point of Lovelock theory, cf. Zanelli (2005). Nevertheless, it is
worth emphasizing that this special curve on the parameter space is not the
one that leads to a single maximally symmetric vacuum in the Einstein-Gauss-
Bonnet theory; in the conventions are using, the maximally symmetric vacua
coincide at $\Lambda\alpha=-5/12$ instead. For generic, fixed values of the
Kähler potential $\mathcal{B}$, we are working with spaces that are not
continuously connected with a maximally symmetric vacuum, regardless of the
precise expression for the $f(r)$; nevertheless, there is still a particular
value of the couplings for which the two branches coincide when both $\sigma$
and $a$ vanish. Since our spaces are not of constant curvature, the
perturbations of the field equations around the solution at this special point
may have support, in contrast to what happens with the perturbations around
the maximally symmetric vacuum when both vacua coincide.

On the curve of the parameter space where $\Lambda\alpha=-1/3$, the radial
coordinate is bounded if $\sigma=1$ according to

$\displaystyle r_{b}\leq
r\leq\left(-\frac{3a^{6}\Lambda}{4}\right)^{\frac{1}{4}}\,,$ (19)

with $\Lambda<0$ and $a\in\mathbb{R}$. When $\sigma=0$, the solution becomes

$\displaystyle f_{\pm}(r)=-\frac{\Lambda
r^{2}}{8}+\frac{2\gamma}{3}\pm\frac{a^{3}}{r^{2}}\sqrt{-\frac{\Lambda}{6}}\,,$
(20)

and it exists only for $\Lambda<0$. The latter implies that both branches have
a bolt at $f_{\pm}(r_{b})=0$ depending on the values of the parameters
$\gamma$ and $a$. In particular, for $\gamma=1$ and $\gamma=0$ we find that
$f_{+}(r)$ and $f_{-}(r)$ have a bolt if $a<0$ and $a>0$, respectively, and
the space is thus completely regular. In contrast, if $\gamma=-1$ with
$\sigma=0$ the metric exhibits a naked curvature singularity at $r=0$.

##  V Topological invariants and renormalized action

In this section, we compute the renormalized Euclidean on-shell action and
Euler characteristic for the Eguchi-Hanson metric (13) with compact base
manifolds $\mathbb{CP}^{2}$, $\mathbb{T}^{4}$, and
$\mathbb{S}^{2}\times\mathbb{S}^{2}$. In the next section we will present the
higher-dimensional generalization of the Eguchi-Hanson metric to arbitrary
even dimension $D$, and, in particular, we will see that the solution in $D=8$
has also non-trivial Pontryagin index, as it happens in $D=4$.

Let us first compute the Euler characteristic $\chi(\mathcal{M}_{6})$ for the
solution presented in Eq. (13). In arbitrary even dimensions $D=2m$, the Euler
theorem states that

$\displaystyle\int_{\mathcal{M}_{2m}}\text{d}^{2m}x\sqrt{|g|}\,\mathcal{E}_{2m}=\left(4\pi\right)^{m}\,m!\,\chi\left(\mathcal{M}_{2m}\right)+\int_{\partial\mathcal{M}_{2m}}\text{d}^{2m-1}x\sqrt{|h|}\,\mathcal{C}_{2m-1}\,,$
(21)

where $\mathcal{E}_{2m}\equiv\mathscr{L}^{(m)}$ [cf. Eq. (7)] is the Euler
density in $D=2m$ dimensions, $\chi\left(\mathcal{M}_{2m}\right)$ is the Euler
characteristic of the $2m$-dimensional manifold $\mathcal{M}_{2m}$, $h$ is the
determinant of the induced metric $h_{\mu\nu}=g_{\mu\nu}-n_{\mu}n_{\nu}$ with
$n^{\mu}$ being a space-like unit vector that is normal to the
($2m-1$)-dimensional boundary $\partial\mathcal{M}_{2m}$, and
$\mathcal{C}_{2m-1}$ is the Chern form on $\partial\mathcal{M}_{2m}$ defined
through the parametric integral as

$\displaystyle\mathcal{C}_{2m-1}$
$\displaystyle=2m\int_{0}^{1}\text{d}s\,\delta^{\mu_{1}\ldots\mu_{2m-1}}_{\nu_{1}\ldots\nu_{2m-1}}K^{\nu_{1}}_{\mu_{1}}\left(\frac{1}{2}\mathcal{R}^{\nu_{2}\nu_{3}}_{\
\
\mu_{2}\mu_{3}}-s^{2}\,K^{\nu_{2}}_{\mu_{2}}K^{\nu_{3}}_{\mu_{3}}\right)\times$
$\displaystyle\quad\ldots\times\left(\frac{1}{2}\mathcal{R}^{\nu_{2m-2}\nu_{2m-1}}_{\
\
\mu_{2m-2}\mu_{2m-1}}-s^{2}\,K^{\nu_{2m-2}}_{\mu_{2m-1}}K^{\nu_{2m-1}}_{\mu_{2m-1}}\right)\,,$
(22)

with $\mathcal{R}^{\mu\nu}_{\ \lambda\rho}$ and
$K_{\mu\nu}=h^{\lambda}_{\mu}\nabla_{\lambda}n_{\nu}$ being the intrinsic and
extrinsic curvature, respectively. These are related with the Riemannian
curvature through the Gauss-Codazzi relation

$\displaystyle\mathcal{R}^{\mu\nu}_{\
\lambda\rho}=h^{\mu}_{\alpha}h^{\nu}_{\beta}h^{\gamma}_{\lambda}h^{\delta}_{\rho}R^{\alpha\beta}_{\
\gamma\delta}+K^{\mu}_{\lambda}K^{\nu}_{\rho}-K^{\mu}_{\rho}K^{\nu}_{\lambda}\,.$
(23)

Focusing on the compact cases of the solution (13), a direct evaluation of the
Euler characteristic in Eq. (21) yields $\chi(\mathcal{M}_{6})=3$,
$\chi(\mathcal{M}_{6})=0$, and $\chi(\mathcal{M}_{6})=4$ for the fibered base
manifolds $\mathbb{CP}^{2}$, $\mathbb{T}^{4}$, and
$\mathbb{S}^{2}\times\mathbb{S}^{2}$, respectively. The latter can be
equivalently obtained by means of the Künneth theorem, since
$\chi(\mathbb{S}^{2}\times\mathbb{S}^{2})=\chi(\mathbb{S}^{2})\times\chi(\mathbb{S}^{2})=2\times
2=4$.

Here, we focus on the curve of the parameter space when $\Lambda\alpha=-1/3$
for the sake of simplicity. For conformally flat boundaries, the renormalized
Euclidean on-shell action can be obtained by appropriately fixing the coupling
constant $\zeta$ of the six-dimensional Euler density [cf. Eq. (9)]. The
latter plays the role of topological counterterms in even-dimensional AlAdS
solutions. Although the Eguchi-Hanson solution in Einstein-Gauss-Bonnet
gravity is not AlAdS, there is still a non-trivial value of the Euler density
coupling that yields a finite result; this corresponds to
$\zeta\Lambda^{2}=1/9$ in Eq. (9). On the other hand, when $\sigma=1$, the
boundary is not conformally flat; nevertheless, the radial coordinate is
naturally bounded according to (19), and so there is no need to include
neither topological counterterms nor conformal completion to render the action
finite. This means that in that case we can set $\zeta=0$ in a natural way.
This is different to what happen when $\sigma=0$, for which the condition
$\zeta\Lambda^{2}=1/9$ is needed for finiteness. Then, considering $\Lambda<0$
and the solutions with compact base manifolds, we obtain

$\displaystyle I_{\rm EGB}$ $\displaystyle=\left\\{\begin{matrix}\frac{4\kappa
a^{9}\beta_{\tau}\pi^{2}}{3r_{b}^{6}}\sqrt{-\frac{6}{\Lambda}}&\mbox{for}&\mathbb{CP}^{2}\,,\\\
\frac{2\kappa
a^{9}\beta_{\tau}\beta_{\theta_{1}}\beta_{\phi_{1}}\beta_{\theta_{2}}\beta_{\phi_{2}}}{27r_{b}^{6}}\sqrt{-\frac{6}{\Lambda}}&\mbox{for}&\mathbb{T}^{4}\,,\\\
-\frac{2\kappa\beta_{\tau}\pi^{2}r_{b}^{2}\left[2r_{b}^{4}+\Lambda
a^{6}-\frac{r_{b}^{4}\Lambda}{12}\sqrt{-(8r_{b}^{4}+6\Lambda
a^{6})}\right]}{\sqrt{-(8r_{b}^{4}+6\Lambda
a^{6})}}&\mbox{for}&\mathbb{S}^{2}\times\mathbb{S}^{2}\,.\end{matrix}\right.$
(24)

The renormalized Euclidean on-shell action can be used in the path integral
approach to study the contribution of these instantons, for instance resorting
to the saddle point approximation Gibbons and Hawking (1977). In the case of
AlAdS spaces, the Euclidean path integral approach permits to explore the
thermodynamics; for instance, it enables to study the well-known phase
transitions that black holes develop in that sector at certain critical
temperature Hawking and Page (1983). A similar behavior may be found for
solutions whose Euclidean version corresponds to AlAdS gravitational
instantons Johnson (2014a, b). Thermodynamics and phase transitions of AlAdS
solutions have also been studied in Lovelock theory, for instance in Cvetic
_et al._ (2002); Nojiri and Odintsov (2002); Cai and Guo (2004); Clunan  _et
al._ (2004); Dehghani and Shamirzaie (2005); Camanho  _et al._ (2014);
Frassino  _et al._ (2014); Aránguiz  _et al._ (2016); Su  _et al._ (2020). In
the case of Eguchi-Hanson solutions, being intrinsically Euclidean and not
AlAdS, the interpretation of the path integral results would be different.
However, it would still be interesting to study whether any sort of phase
transitions occurs in this sector. Instantons can contribute to the vacuum
persistence amplitude in gauge theory Osborn (1981), and a similar computation
can be performed in the case of gravitational instantons, cf. Strominger  _et
al._ (1984). In the latter work, the authors provided a general formula for
the one-loop determinant in terms of the renormalization group invariant
masses, the volume of space, and topological invariants, by performing linear
perturbations around the instantonic background. We believe these are
interesting applications that can be studied using the analytic solutions
found here.

In the next section, we generalize the Eguchi-Hanson instanton of Einstein-
Gauss-Bonnet gravity to arbitrary even dimensions. Moreover, we show how the
solution is modified in presence of cubic curvature terms in the Lovelock
series.

##  VI Generalization to higher dimensions

Let us extend our solution to higher dimension. The general even-dimensional
case, say $D=2m$ with $m\geq 3$, can be solved analytically by considering
different topologies of the base manifold. This was done in Sec. IV for six
dimensions as a particular case, and here we consider the metric Ansatz (1)
with base manifold of topology $\mathbb{CP}^{k}$, $\mathbb{CH}^{k}$,
$(\mathbb{T}^{2})^{k}$, $(\mathbb{S}^{2})^{k}$, and $(\mathbb{H}^{2})^{k}$;
the last three are the product of $k=m-1$ constant curvature 2-manifolds.
Compact base manifolds with negative constant curvature can also be
considered, for example, by taking spaces locally equivalent to
$(\mathbb{H}^{2})^{k}$ constructed out of quotients of the hyperbolic space by
discrete subgroups. The Kähler potential and line element for each of these
$(2m-2)$-dimensional base manifolds is given below. The field equations of
Einstein-Gauss-Bonnet gravity (11) are solved analytically by the line element
(1) with the metric function

$\displaystyle f_{\pm}(r)$
$\displaystyle=\frac{2}{m}\Bigg{\\{}\gamma+\frac{r^{2}}{16\alpha(m-2)}\left[1\pm\sqrt{1+8(m-2)\left(\frac{\Lambda\alpha
m}{m^{2}-1}+\frac{4\alpha
a^{2m}}{r^{2(m+1)}}-\frac{32\sigma\alpha^{2}(m-2)}{r^{4}(m-1)}\right)}\,\right]\Bigg{\\}}\,.$
(25)

Here $a$ is an integration constant, while $\gamma$ and $\sigma$ are
parameters that depend on the topology and they admit the same geometric
interpretation as the one given above Table 14. In particular, for the
topologies considered in this Section, their values are

| $\;\;\;\mathbb{CP}^{k}\;\;\;$ | $\;\;\;\mathbb{CH}^{k}\;\;\;$ | $\;(\mathbb{T}^{2})^{k}\;$ | $\;(\mathbb{S}^{2})^{k}\;$ | $\;(\mathbb{H}^{2})^{k}\;$
---|---|---|---|---|---
$\;\;\;\gamma\;\;\;$ | $1$ | $-1$ | $0$ | $1$ | $-1$
$\;\;\;\sigma\;\;\;$ | $0$ | $0$ | $0$ | $1$ | $1$
(26)

Kähler spaces naturally allow for circle fiber bundles to be constructed on
them. This is relevant for the construction of the higher-dimensional versions
of the Eguchi-Hanson space discussed here. Kähler manifolds are described by a
Riemannian metric with associated line element $\text{d}\Sigma^{2}$ and a
compatible symplectic closed 2-form $\Omega$. It is convenient to parameterize
the symplectic form by a potential 1-form $\cal B$ defined by
$\Omega=\text{d}\cal B$. For our purposes, we require the Kähler base
manifolds to be, in addition, Einstein spaces. In two dimensions all spaces
trivially fulfill these conditions; in higher dimensions, however, this is
actually a stringent constraint. One way of constructing higher-dimensional
Kähler-Einstein manifolds in arbitrary even dimension $2k$ is as the direct
product of $k$ 2-dimensional spaces of equal curvature. For example, the
product spaces $(\mathbb{T}^{2})^{k}$, $(\mathbb{S}^{2})^{k}$ and
$(\mathbb{H}^{2})^{k}$ meet these requirements whenever their geometries are
of the form detailed in the following Table:

Space |  ${\cal B}_{i}$ |  $\text{d}\Sigma_{i}^{2}$ |  $\theta_{i}$ range |  $\phi_{i}$ range
---|---|---|---|---
$(\mathbb{T}^{2})^{k}$ |  $\theta_{i}\text{d}\phi_{i}$ |  $\text{d}\theta_{i}^{2}+\text{d}\phi_{i}^{2}$ |  $[0,\beta_{\theta_{i}}]$ |  $[0,\beta_{\phi_{i}}]$
$(\mathbb{S}^{2})^{k}$ |  $\cos\theta_{i}\text{d}\phi_{i}$ |  $\hskip 56.9055pt\text{d}\theta_{i}^{2}+\sin^{2}\theta_{i}\text{d}\phi_{i}^{2}$ |  $[0,\pi]$ |  $[0,2\pi]$
$(\mathbb{H}^{2})^{k}$ |  $\cosh\theta_{i}\text{d}\phi_{i}$ |  $\hskip 56.9055pt\text{d}\theta_{i}^{2}+\sinh^{2}\theta_{i}\text{d}\phi_{i}^{2}$ |  $(-\infty,\infty)$ |  $[0,2\pi]$

This Table shows a list of Kähler-Einstein spaces with symplectic potential
${\cal B}=\sum^{k}_{i=1}{\cal B}_{i}$ and a compatible Riemannian metric
$\text{d}\Sigma^{2}=\sum^{k}_{i=1}\text{d}\Sigma^{2}_{i}$. For $k>1$, even-
dimensional hyperspheres are Einstein spaces but not Kähler manifolds. In
contrast, complex projective spaces $\mathbb{CP}^{k}$ are Kähler-Einstein
manifolds for arbitrary $k$ when they are equipped with the Fubini-Study
metric. Complex hyperbolic spaces with the Bergman metric also fulfill the
geometric desiderata for our construction. Part of our calculations is
benefited from considering the iterative construction of complex projective
spaces Hoxha  _et al._ (2000). In this case, both the metric and the potential
are written recursively; namely

$\displaystyle\mathcal{B}_{(k)}$
$\displaystyle=(k+1)\sin^{2}\psi_{k}\left(\text{d}\phi_{k}+\frac{1}{k}\mathcal{B}_{(k-1)}\right),$
(27a) $\displaystyle\text{d}\Sigma_{(k)}^{2}$
$\displaystyle=2(k+1)\bigg{[}\text{d}\psi_{k}^{2}+\sin^{2}\psi_{k}\cos^{2}\psi_{k}\left(\text{d}\phi_{k}+\frac{1}{k}\mathcal{B}_{(k-1)}\right)^{2}+\frac{1}{2k}\sin^{2}\psi_{k}\text{d}\Sigma_{(k-1)}^{2}\bigg{]}\,,$
(27b)

where $0\leq\psi_{k}\leq\pi/2$ and $0\leq\phi_{k}\leq 2\pi$. The base case
($k=1$) is provided by the sphere. The iterative formulae above take advantage
of the fact that $\mathbb{CP}^{k}$ admits a foliation by $\mathbb{S}^{2k-1}$,
modulo two points. The complex hyperbolic space $\mathbb{CH}^{k}$ admits a
similar slicing; namely,

$\displaystyle\mathcal{B}$
$\displaystyle=(k+1)\sinh^{2}\psi_{k}\left(\text{d}\phi_{k}+\frac{1}{k}\mathcal{B}_{(k-1)}\right),$
(28a) $\displaystyle\text{d}\Sigma^{2}$
$\displaystyle=2(k+1)\bigg{[}\text{d}\psi_{k}^{2}+\sinh^{2}\psi_{k}\cosh^{2}\psi_{k}\left(\text{d}\phi_{k}+\frac{1}{k}\mathcal{B}_{(k-1)}\right)^{2}+\frac{1}{2k}\sinh^{2}\psi_{k}\text{d}\Sigma_{(k-1)}^{2}\bigg{]}\,,$
(28b)

where $-\infty<\psi_{k}<\infty$ but all else is as in Eq. (27).

As in the lower-dimensional examples, the bolt of the $2m$-dimensional
solution we are constructing is defined as the ($2m-2$)-dimensional set of
fixed points, located at the locus $f_{\pm}(r_{b})=0$. This condition,
alongside the absence of conical singularities, leads to a space that is
geodesically complete. The branch $f_{-}(r)$ is continuously connected to the
Einstein solution in arbitrary even dimensions [cf. Eq. (3)] in the limit
$\alpha\to 0$. Additionally, defining
$\Delta\equiv(m^{2}-1)\left[8\Lambda\alpha m(m-2)+m^{2}-1\right]$, we find
that the asymptotic behavior of this metric is given by

$\displaystyle f(r)$
$\displaystyle\simeq\frac{\left[1\pm\sqrt{\frac{\Delta}{(m^{2}-1)^{2}}}\right]r^{2}}{8\alpha
m(m-2)}+\frac{2\gamma}{m}\mp\frac{16\alpha\sigma
m(m+1)}{(m-2)\sqrt{\Delta}\,r^{2}}\mp\frac{1024\alpha^{3}\sigma^{2}m^{3}(m+1)^{2}}{(m-2)\sqrt{\Delta}\,r^{6}}+\mathcal{O}(r^{-10})\,,$

as $r\to\infty$. The asymptotia of the branch $f_{-}(r)$ matches the Einstein
gravity solution at large $r$ for any Kähler-Einstein base manifold. Solution
(25) is the generalized Eguchi-Hanson space in arbitrary even dimension
$D=2m$, and in presence of quadratic curvature terms.

##  VII Generalizations with cubic terms and higher

The solution we have just presented can be extended to Lovelock theory with
higher-curvature terms, cf. Wheeler (1986a, b). As an example, in eight
dimensions ($m=4$) consider a 6-torus base manifold $\mathbb{T}^{6}$ (i.e.
$k=3$). For this ansatz, we can solve the equations of motion for the
Lagrangian (7) with $p=0,1,2,3$ explicitly. That is to say, consider, in
addition to the quadratic terms in (9)-(10), the following cubic term in the
action

$\displaystyle I^{(3)}[g]=$
$\displaystyle\frac{\alpha_{3}\kappa}{8}\int_{\mathcal{M}}\text{d}^{2m}x\sqrt{|g|}\,\left(R^{3}-12RR_{\mu\nu}R^{\mu\nu}+16R_{\mu\nu}R^{\mu}_{\,\,\rho}R^{\nu\rho}\right.$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
+24R_{\mu\nu}R_{\rho\eta}R^{\mu\nu\rho\eta}+3RR_{\mu\nu\rho\eta}R^{\mu\nu\rho\eta}-24R_{\mu\nu}R^{\mu}_{\,\rho\eta\sigma}R^{\nu\rho\eta\sigma}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
\left.+4R_{\mu\nu\rho\eta}R^{\mu\nu\alpha\beta}R^{\rho\eta}_{\,\,\,\,\alpha\beta}-8R_{\mu\nu\rho}^{\,\,\,\,\,\,\,\,\,\eta}R^{\mu\alpha\rho\beta}R^{\nu}_{\,\,\alpha\eta\beta}\right)\,.$

It can easily be checked that a polynomial equation implicitly determines the
metric function $f(r)$. More precisely, in $D=8$, the field equations now take
the form $G_{\mu\nu}+\Lambda
g_{\mu\nu}+\alpha_{2}H_{\mu\nu}+\alpha_{3}M_{\mu\nu}=0$, where $H_{\mu\nu}$ is
defined in Eq. (12) and $M_{\mu\nu}$ comes from the variation of (VII); namely

$\displaystyle M^{\mu}_{\nu}$
$\displaystyle=-\frac{1}{16}\delta^{\mu\lambda_{2}\ldots\lambda_{7}}_{\nu\rho_{2}\ldots\rho_{7}}R_{\
\ \lambda_{2}\lambda_{3}}^{\rho_{2}\rho_{3}}R_{\ \
\lambda_{4}\lambda_{5}}^{\rho_{4}\rho_{5}}R_{\ \
\lambda_{6}\lambda_{7}}^{\rho_{6}\rho_{7}}\,.$ (30)

This yields the Wheeler type polynomial Wheeler (1986a, b)

$\displaystyle
384\alpha_{3}r^{4}f(r)^{3}-32\alpha_{2}f(r)^{2}r^{6}+r^{8}f(r)+\frac{\Lambda}{30}r^{10}+a^{8}=0\,,$
(31)

where $a$ is an integration constant. To see how this polynomial compares to
others like it in the literature, recall that the Eguchi-Hanson metric closely
resembles that of Taub-NUT, and this similarity carries over to their
respective generalizations. However, in spite of these parallels Eq. (31) is
quite different from those found for Taub-NUT, see Ref. Corral  _et al._
(2019).

Furthermore, Eq. (31) is a polynomial equation that can be solved explicitly;
however, its generic solution is cumbersome and not particularly illuminating.
In order to have a clearer picture of the cubic solution, let us consider a
particular case that enables us to visualize the form of the metric: Consider
$\alpha_{3}={8\alpha_{2}^{2}}/{9}$, for which we find

$\displaystyle
f(r)=\frac{r^{2}}{32\alpha_{2}}-\frac{r^{2}}{32\alpha_{2}}\left(1+\frac{16\alpha_{2}\Lambda}{5}+\frac{96\alpha_{2}a^{8}}{r^{10}}\right)^{\frac{1}{3}}\,.$
(32)

This is an exact 8-dimensional solution of the cubic theory that behaves as
the solution of $8$-dimensional Einstein gravity in the small $\alpha_{2}$
limit; namely $f(r)\simeq-\frac{\Lambda
r^{2}}{30}-\frac{a^{8}}{r^{8}}+\mathcal{O}(\alpha_{2})$. It is asymptotically
equivalent to (1) with $\text{d}\Sigma^{2}$ being the flat metric on
$\mathbb{T}^{6}$. Function (32) has its root(s) at
$r_{b}=(-30a^{8}/\Lambda)^{1/10}$. This is a concrete example of exact,
analytic solution of the cubic Lovelock theory that asymptotically approaches
a higher-dimensional generalization of Eguchi-Hanson space in presence of
cosmological constant. Further examples can easily be obtained by following
the method we described in this paper, i.e. by solving a polynomial of degree
$n$, analogous to (31), when the $2m$-dimensional extension of the
$n$-dimensional Euler characteristics of order $\mathcal{O}(R^{n})$, with
$n\leq m-1$, are added to the gravity action.

###### Acknowledgements.

We thank to Giorgos Anastasiou, Ignacio J. Araya, Olivera Mišković, and
Rodrigo Olea for insightful comments and discussions. The work of C. C. is
partially supported by Agencia Nacional de Investigación y Desarrollo (ANID)
through FONDECYT grants No 11200025 and 1210500. D. F. is supported by a
CONACYT postdoctoral fellowship. This work has been partially funded by
CONACYT Grant No. A1-S-11548 and by Dirección de Postgrado (UdeC), through the
grant UCO1866. G. G. is supported by CONICET and ANPCyT grants PIP-1109-2017,
PICT-2019-00303. J.O. is partially supported by FONDECYT grants 1221504 and
1210635 and by Proyecto de cooperación internacional 2019/13231-7 FAPESP/ANID.

## References

  * Eguchi and Hanson (1978) T. Eguchi and A. J. Hanson, Phys. Lett. B 74, 249 (1978).
  * Eguchi and Hanson (1979a) T. Eguchi and A. J. Hanson, Annals Phys. 120, 82 (1979a).
  * Calabi (1979) E. Calabi, Annales Scientifiques de l’École Normale Supérieure 12, 269 (1979).
  * Hawking (1977) S. W. Hawking, Phys. Lett. A 60, 81 (1977).
  * Eguchi and Freund (1976) T. Eguchi and P. G. O. Freund, Phys. Rev. Lett. 37, 1251 (1976).
  * Gibbons and Hawking (1978) G. W. Gibbons and S. W. Hawking, Phys. Lett. B 78, 430 (1978).
  * Page (1978) D. N. Page, Phys. Lett. B 78, 249 (1978).
  * Hawking and Pope (1978) S. W. Hawking and C. N. Pope, Nucl. Phys. B 146, 381 (1978).
  * Gibbons and Hawking (1979) G. W. Gibbons and S. W. Hawking, Commun. Math. Phys. 66, 291 (1979).
  * Eguchi and Hanson (1979b) T. Eguchi and A. J. Hanson, Gen. Rel. Grav. 11, 315 (1979b).
  * Eguchi  _et al._ (1980) T. Eguchi, P. B. Gilkey, and A. J. Hanson, Phys. Rept. 66, 213 (1980).
  * Pedersen (1985) H. Pedersen, Class. and Quan. Gravity 2, 579 (1985).
  * Higashijima  _et al._ (2002a) K. Higashijima, T. Kimura, and M. Nitta, Nucl. Phys. B 623, 133 (2002a).
  * Higashijima  _et al._ (2002b) K. Higashijima, T. Kimura, and M. Nitta, Annals Phys. 296, 347 (2002b).
  * Higashijima  _et al._ (2002c) K. Higashijima, T. Kimura, and M. Nitta, Nucl. Phys. B 645, 438 (2002c).
  * Clarkson and Mann (2006a) R. Clarkson and R. B. Mann, Class. Quant. Grav. 23, 1507 (2006a).
  * Clarkson and Mann (2006b) R. Clarkson and R. B. Mann, Phys. Rev. Lett. 96, 051104 (2006b).
  * Giribet and Santillan (2007) G. Giribet and O. P. Santillan, Commun. Math. Phys. 275, 373 (2007).
  * Wong and Mann (2012) A. W. C. Wong and R. B. Mann, Phys. Rev. D 85, 046010 (2012).
  * Lovelock (1971) D. Lovelock, J. Math. Phys. 12, 498 (1971).
  * Lanczos (1938) C. Lanczos, Annals Math. 39, 842 (1938).
  * Lovelock (1972) D. Lovelock, J. Math. Phys. 13, 874 (1972).
  * Hendi and Dehghani (2008) S. H. Hendi and M. H. Dehghani, Phys. Lett. B666, 116 (2008).
  * Hendi  _et al._ (2012) S. H. Hendi, R. B. Mann, N. Riazi, and B. Eslam Panah, Phys. Rev. D 86, 104034 (2012).
  * Bueno  _et al._ (2018) P. Bueno, P. A. Cano, R. A. Hennigar, and R. B. Mann, JHEP 10, 095 (2018).
  * Corral  _et al._ (2019) C. Corral, D. Flores-Alfonso, and H. Quevedo, Phys. Rev. D100, 064051 (2019).
  * Corral  _et al._ (2021) C. Corral, G. Giribet, and R. Olea, Phys. Rev. D 104, 064026 (2021).
  * Page and Pope (1987) D. N. Page and C. N. Pope, Class. Quant. Grav. 4, 213 (1987).
  * Belavin  _et al._ (1975) A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Y. S. Tyupkin, Phys. Lett. B 59, 85 (1975).
  * Misner (1963) C. W. Misner, J. Math. Phys. 4, 924 (1963).
  * Hawking and Hunter (1999) S. W. Hawking and C. J. Hunter, Phys. Rev. D 59, 044025 (1999).
  * Chen and Zhang (2021) J. Chen and X. Zhang, J. Geom. Phys. 161, 104010 (2021).
  * Atiyah  _et al._ (1975) M. F. Atiyah, V. K. Patodi, and I. M. Singer, Math. Proc. Cambridge Phil. Soc. 77, 43 (1975).
  * Atiyah  _et al._ (1976a) M. F. Atiyah, V. K. Patodi, and I. M. Singer, Math. Proc. Cambridge Phil. Soc. 79, 71 (1976a).
  * Atiyah  _et al._ (1976b) M. F. Atiyah, V. K. Patodi, and I. M. Singer, Math. Proc. Cambridge Phil. Soc. 78, 405 (1976b).
  * Henneaux  _et al._ (1990) M. Henneaux, C. Teitelboim, and J. Zanelli, Nucl. Phys. B 332, 169 (1990).
  * Chern (1945) S.-s. Chern, The Annals of Mathematics 46, 674 (1945).
  * Aros  _et al._ (2000a) R. Aros, M. Contreras, R. Olea, R. Troncoso, and J. Zanelli, Phys. Rev. Lett. 84, 1647 (2000a).
  * Aros  _et al._ (2000b) R. Aros, M. Contreras, R. Olea, R. Troncoso, and J. Zanelli, Phys. Rev. D 62, 044002 (2000b).
  * Mora  _et al._ (2004) P. Mora, R. Olea, R. Troncoso, and J. Zanelli, JHEP 06, 036 (2004).
  * Kofinas and Olea (2006) G. Kofinas and R. Olea, Phys. Rev. D 74, 084035 (2006).
  * Olea (2007) R. Olea, JHEP 04, 073 (2007).
  * Arenas-Henriquez  _et al._ (2019) G. Arenas-Henriquez, R. B. Mann, O. Miskovic, and R. Olea, Phys. Rev. D 100, 064038 (2019).
  * Ciambelli  _et al._ (2021) L. Ciambelli, C. Corral, J. Figueroa, G. Giribet, and R. Olea, Phys. Rev. D 103, 024052 (2021).
  * Anastasiou  _et al._ (2020) G. Anastasiou, O. Miskovic, R. Olea, and I. Papadimitriou, JHEP 08, 061 (2020).
  * Anastasiou  _et al._ (2021a) G. Anastasiou, I. J. Araya, and R. Olea, JHEP 01, 134 (2021a).
  * Anastasiou  _et al._ (2021b) G. Anastasiou, I. J. Araya, C. Corral, and R. Olea, JHEP 2021, 156 (2021b).
  * Miskovic and Olea (2009) O. Miskovic and R. Olea, Phys. Rev. D 79, 124020 (2009).
  * Araneda  _et al._ (2016) R. Araneda, R. Aros, O. Miskovic, and R. Olea, Phys. Rev. D 93, 084022 (2016).
  * Boulware and Deser (1985) D. G. Boulware and S. Deser, Phys. Rev. Lett. 55, 2656 (1985).
  * Cai (2002) R.-G. Cai, Phys. Rev. D 65, 084014 (2002).
  * Dotti  _et al._ (2007a) G. Dotti, J. Oliva, and R. Troncoso, Phys. Rev. D 75, 024002 (2007a).
  * Dotti  _et al._ (2007b) G. Dotti, J. Oliva, and R. Troncoso, Phys. Rev. D 76, 064038 (2007b).
  * Dotti and Gleiser (2005) G. Dotti and R. J. Gleiser, Phys. Lett. B 627, 174 (2005).
  * Zanelli (2005) J. Zanelli, in  _7th Mexican Workshop on Particles and Fields_ (2005) arXiv:hep-th/0502193 .
  * Gibbons and Hawking (1977) G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2752 (1977).
  * Hawking and Page (1983) S. W. Hawking and D. N. Page, Commun. Math. Phys. 87, 577 (1983).
  * Johnson (2014a) C. V. Johnson, Class. Quant. Grav. 31, 235003 (2014a).
  * Johnson (2014b) C. V. Johnson, Class. Quant. Grav. 31, 225005 (2014b).
  * Cvetic  _et al._ (2002) M. Cvetic, S. Nojiri, and S. D. Odintsov, Nucl. Phys. B 628, 295 (2002).
  * Nojiri and Odintsov (2002) S. Nojiri and S. D. Odintsov, Phys. Rev. D 66, 044012 (2002).
  * Cai and Guo (2004) R.-G. Cai and Q. Guo, Phys. Rev. D 69, 104025 (2004).
  * Clunan  _et al._ (2004) T. Clunan, S. F. Ross, and D. J. Smith, Class. Quant. Grav. 21, 3447 (2004).
  * Dehghani and Shamirzaie (2005) M. H. Dehghani and M. Shamirzaie, Phys. Rev. D 72, 124015 (2005).
  * Camanho  _et al._ (2014) X. O. Camanho, J. D. Edelstein, G. Giribet, and A. Gomberoff, Phys. Rev. D 90, 064028 (2014).
  * Frassino  _et al._ (2014) A. M. Frassino, D. Kubiznak, R. B. Mann, and F. Simovic, JHEP 09, 080 (2014).
  * Aránguiz  _et al._ (2016) L. Aránguiz, X.-M. Kuang, and O. Miskovic, Phys. Rev. D 93, 064039 (2016).
  * Su  _et al._ (2020) B.-Y. Su, Y.-Y. Wang, and N. Li, Eur. Phys. J. C 80, 305 (2020).
  * Osborn (1981) H. Osborn, Annals Phys. 135, 373 (1981).
  * Strominger  _et al._ (1984) A. Strominger, G. T. Horowitz, and M. J. Perry, Nucl. Phys. B 238, 653 (1984).
  * Hoxha  _et al._ (2000) P. Hoxha, R. R. Martinez-Acosta, and C. N. Pope, Class. Quant. Grav. 17, 4207 (2000).
  * Wheeler (1986a) J. T. Wheeler, Nucl. Phys. B 268, 737 (1986a).
  * Wheeler (1986b) J. T. Wheeler, Nucl. Phys. B 273, 732 (1986b).


††thanks: These authors contributed equally††thanks: These authors contributed
equally

# Supplemental Information for “Transient amplification of broken
symmetry in elastic snap-through”

Qiong Wang  Department of Mechanical Science and Engineering, University of
Illinois at Urbana-Champaign, Urbana, IL 61801  Andrea Giudici  Mathematical
Institute, University of Oxford, Woodstock Rd, Oxford, OX2 6GG, UK  Weicheng
Huang  School of Mechanical Engineering, Southeast University, Nanjing 211189,
China  Yuzhe Wang  Singapore Institute of Manufacturing Technology, Agency for
Science, Technology and Research, Singapore, 138634, Singapore  Mingchao Liu
Mathematical Institute, University of Oxford, Woodstock Rd, Oxford, OX2 6GG,
UK  School of Mechanical and Aerospace Engineering, Nanyang Technological
University, Singapore 639798, Singapore  Department of Mechanical Engineering,
University of Birmingham, Birmingham, B15 2TT, UK  Sameh Tawfick  Department
of Mechanical Science and Engineering, University of Illinois at Urbana-
Champaign, Urbana, IL 61801  Dominic Vella<EMAIL_ADDRESS>Mathematical Institute, University of Oxford, Woodstock Rd, Oxford, OX2 6GG,
UK

In this Supplementary Information, we present further details of the
experimental setup (§I), the analysis that leads to the bifurcation diagram in
fig. 1c of the main paper (§II) and the numerical simulations (§III). The bulk
of the SI, however, focuses on the development of the toy model (§IV),
including giving further details of the various constants and the
simplifications that arise when the angle $\alpha\ll 1$.

##  I Experimental setup

During the experiment, a strip of spring steel (from McMASTER-CARR, Wear-
Resistant 1095 Spring Steel strip, modulus $E\approx 200\mathrm{~{}GPa}$,
thickness $b=0.05\mathrm{~{}mm}$, width $w=25\mathrm{~{}mm}$) of length
$L=108\mathrm{~{}mm}$ was formed into an arch by clamping its ends
symmetrically at an angle $\alpha=\pi/12\mathrm{~{}rad}$ to the horizontal
with two holders.

The loading setup comprises a linear guide, a linear motorized stage, and a
stopper as shown in fig. S1a. This ensures that the holders can be moved
relative to one another: the holder on the right side is fixed in the
laboratory frame, while the holder on the left side is connected to the
motorized linear stage via a piece of double-sided tape. As the stage moves,
the holder moves with it until it hits the stopper, at which point it detaches
from the moving stage. (The stopper is placed just beyond the position at
which snap-through is expected. This position is determined by moving the left
holder quasistatically until snap-through is observed.) Videos of the snap-
through process were taken with a high-speed camera (Fastec Hispec 5). The
images in the video contained $1184\mathrm{~{}pixels}\times
424\mathrm{~{}pixels}$; the resolution is $96\mathrm{~{}pix/in}\times
96\mathrm{~{}pix/in}$.

The driving frequency and acceleration of the motor that drives the stage are
pre-programmed; the frequency ranged from $10\mathrm{~{}Hz}$ to
$10^{5}\mathrm{~{}Hz}$ while the acceleration time was fixed at the minimum
setting ($0.01\mathrm{~{}s}$). Within one test, the motor first accelerates,
taking the left holder with it and decreasing the end shortening $\Delta L$.
As $\Delta L$ further decreases and passes the critical value where the snap-
through happens, the holder hits the stopper and detaches from the motor while
the motor continues to move before decelerating.

All videos are post-processed in Matlab. The grayscale images are binarized
and the strip is identified using the bwboundaries function. After extracting
the shape of the strip from the videos, the shape data of the strip is
smoothed by Lowess Smoothing [1]. With the shape of the strip at each time
stamp determined, we evaluate the geometrical parameters for analysis by
calculating the slope of the arch at specific locations. In particular, the
inclination angle at the centre, $\theta(t)$, is approximated by determining
the shape of the strip within the 50 pixels on either side of the center-line.
The amplification of the inclination angle $\theta$ is determined by the ratio
between the amplitude of oscillations in $\theta$ before and after snap-
through. The amplitude of oscillations before snap-through is determined from
the standard deviation $\sigma$ to be $\sqrt{2}\sigma$. The amplitude after
the snap-through is given by the peak value of the inclination angle $\theta$
after subtracting the mean value of the inclination angle before the snap-
through.

Figure S1: (a) Top view of the experimental setup showing the steel strip
(lower right). (b) The geometry of a strip with clamped-clamped boundary
conditions. (c) Bifurcation diagram showing the height of the arch center,
$h=\bar{\zeta}(0)$, as a function of $\mu$ in equilibrium. Symmetric modes
with the lowest three compressive forces are plotted as the black (natural
mode), brown (inverted mode) and grey curves. The asymmetric mode with the
lowest compressive force is shown as dashed cyan line. Stable solutions are
shown as full curves, unstable as dashed.

##  II Bifurcation behavior

In this section, we describe in more detail the analysis that leads to the
bifurcation diagram shown in fig. 1c of the main text.

###  II.1 Governing equations and non-dimensionalization

Consider an initially straight, thin strip of length $L$, thickness $h$ and
width $b$ clamped at its two ends with angle $\alpha$ as shown in fig. S1b).
We assume that the material is isotropic and homogeneous with Young’s modulus
$E$ and density $\rho$. Since the strip is slender ($h,b\ll l$), we henceforth
assume that it is unshearable, inextensible and ignore any damping; in this
limit, small transverse deflections, $\bar{\zeta}(\bar{x},\bar{t})$, must
satisfy the dynamic beam equation

$\displaystyle\rho
bh\frac{\partial^{2}\bar{\zeta}}{\partial\bar{t}^{2}}+B\frac{\partial^{4}\bar{\zeta}}{\partial\bar{x}^{4}}+T\frac{\partial^{2}\bar{\zeta}}{\partial\bar{x}^{2}}=0$
(S1)
$\displaystyle\int^{\frac{L}{2}}_{-\frac{L}{2}}\left(\frac{\partial\bar{\zeta}}{\partial\bar{x}}\right)^{2}~{}d\bar{x}=2\frac{\Delta
L}{L}$ (S2)
$\displaystyle\bar{\zeta}\left(\pm\frac{L}{2},\bar{t}\right)=0,\quad\bar{\zeta}_{x}\left(\pm\frac{L}{2},\bar{t}\right)=\mp\alpha,$
(S3)

where $\bar{\zeta}$ is the vertical displacement of the strip, $\alpha$ is the
clamping angle, $T$ is then compressive force applied to the strip and
$B=\frac{Eh^{3}b}{12}$ is the bending stiffness. The vertical displacement at
the center of the strip is denoted by $h$ in the main text.

To non-dimensionalize the equations, we introduce the following dimensionless
variables

$\displaystyle\zeta=\frac{\bar{\zeta}}{L},\quad
x=\frac{\bar{x}}{L},\quad\tau^{2}=\frac{TL^{2}}{B},\quad
t=\frac{\bar{t}}{t^{*}},\quad t^{*}=L^{2}\sqrt{\frac{\rho bh}{B}}$ (S4)

where we use $\tau^{2}$ to denote the dimensionless compressive force for
later notational convenience.

Then eqn. (S1) can be re-written in dimensionless form as

$\displaystyle\frac{\partial^{2}\zeta}{\partial
t^{2}}+\frac{\partial^{4}\zeta}{\partial
x^{4}}+\tau^{2}\frac{\partial^{2}\zeta}{\partial x^{2}}=0$ (S5)
$\displaystyle\int^{\frac{1}{2}}_{-\frac{1}{2}}\left(\frac{\partial\zeta}{\partial
x}\right)^{2}~{}dx=2\frac{\Delta L}{L}$ (S6)
$\displaystyle\zeta\left(\pm\frac{1}{2},t\right)=0,\quad\zeta_{x}\left(\pm\frac{1}{2},t\right)=\mp\alpha.$
(S7)

###  II.2 Static bifurcation behavior

In equilibrium, eqn. (S5)–(S7) can be simplified to

$\displaystyle\frac{d^{4}\zeta}{dx^{4}}+\tau^{2}\frac{d^{2}\zeta}{dx^{2}}=0$
(S8)
$\displaystyle\int^{\frac{1}{2}}_{-\frac{1}{2}}\left(\frac{d\zeta}{dx}\right)^{2}~{}dx=2\frac{\Delta
L}{L}$ (S9)
$\displaystyle\zeta\left(\pm\frac{1}{2}\right)=0,\quad\left.\frac{d\zeta}{dx}\right|_{x=\pm\frac{1}{2}}=\mp\alpha.$
(S10)

The above equations have the general solutions

$\displaystyle\zeta=A\sin{\tau x}+B\cos{\tau x}+C+Dx$ (S11)

where the constants $A,B,C$ and $D$ can be determined from the boundary
conditions eqn. (S10), which now read:

$\displaystyle\begin{split}0&=A\sin{\frac{\tau}{2}}+B\cos{\frac{\tau}{2}}+C+\frac{D}{2}\\\
0&=-A\sin{\frac{\tau}{2}}+B\cos{\frac{\tau}{2}}+C-\frac{D}{2}\\\
\alpha&=A\tau\cos{\frac{\tau}{2}}-B\tau\sin{\frac{\tau}{2}}+D\\\
-\alpha&=A\tau\cos{\frac{\tau}{2}}+B\tau\sin{\frac{\tau}{2}}+D\end{split}$
(S12)

There are two types of solutions to eqn. (S12): symmetric solutions, in which
$A=D=0$, and asymmetric solutions ($A,D\neq 0$), which are only possible if
$\tan{\frac{\tau}{2}}=\frac{\tau}{2}$.

####  II.2.1 Symmetric solutions

For the symmetric case, the solution is

$\displaystyle\zeta_{{\mathrm{eqm}}}^{s}=\frac{\alpha}{\tau}\left(\frac{\cos{\tau
x}}{\sin{\frac{\tau}{2}}}-\cot\frac{\tau}{2}\right)$ (S13)

Here $\tau$ is not known a priori, but must be found from the condition of
inextensibility eqn. (S9), which gives

$\displaystyle\mu=\frac{\Delta
L}{\alpha^{2}L}=\frac{\tau-\sin{\tau}}{4\tau\sin^{2}{\frac{\tau}{2}}}$ (S14)

we shall notice that the value of $\tau$ is solely dependent on the single
parameter $\mu$, which emerges as the ratio between the relative compression
applied and the clamping angle imposed to the strip.

Equations (S13) and (S14) allow the bifurcation diagram to be plotted
parametrically in terms of $\tau$, for example showing $h=\zeta(0)$ as a
function of $\mu$. This is what is plotted in fig. 1c of the main paper.

Fig. S1 depicts $\mu(\tau)$ from (S14) with three lowest solutions of $\tau$
by the black, brown and grey curves. The natural mode corresponds to the
lowest value of $\tau$ (black curve) while the inverted mode corresponds to
the second lowest value of $\tau$ (brown curve). The plot reveals that there
is a local minimum of $\mu$, which is readily calculated to occur at
$\tau\approx 9.203$ and corresponds to $\mu={\mu_{\ast}}\approx 0.247$. For
$\mu>{\mu_{\ast}}\approx 0.247$, both the natural mode and the inverted mode
exists while for $\mu<{\mu_{\ast}}$, only the natural mode exists — a saddle
node bifurcation occurs at $\mu={\mu_{\ast}}$.

####  II.2.2 Asymmetric solutions

Asymmetric solutions exist only if $\tau=\tau_{*}$ such that
$\tau_{*}/2=\tan(\tau_{*}/2)$. We find numerically that $\tau_{*}\approx
8.987$ and the solution may be written

$\displaystyle\zeta_{\mathrm{eqm}}^{a}(x)=\frac{\alpha}{\tau_{*}}\left(\frac{\cos{\tau_{*}x}}{\sin{\frac{\tau_{*}}{2}}}-\cot\frac{\tau_{*}}{2}\right)+\frac{\alpha}{\tau_{*}}D\left(\sin{\tau_{*}x}-2x\sin{\frac{\tau_{*}}{2}}\right).$
(S15)

Note that the symmetric portion of (S15) is precisely that of (S13), albeit
with $\tau$ replaced by $\tau_{*}$. The size of the asymmetry, encoded by $D$,
is determined from the condition of inextensibility eqn. (S9); we find that

$\displaystyle D=\frac{\sqrt{\frac{4\Delta
L}{L\alpha^{2}}-1}}{\sin{\frac{\tau_{*}}{2}}}.$ (S16)

Note that this constant is only real when $\Delta L/(L\alpha^{2})=\mu>1/4$ —
this solution vanishes for $\mu\leq 1/4$.

The bifurcation diagram fig. S1c) is completed by adding the asymmetric mode
(cyan curve), noting that the mid-point displacement
$\zeta(0)=\frac{\alpha}{\tau_{*}}\left(\frac{1-\cos\frac{\tau_{*}}{2}}{\sin\frac{\tau_{*}}{2}}\right)$
is constant, independent of $\mu\geq 1/4$.

###  II.3 Stability analysis

Having determined the various equilibrium solutions that exist at a given
$\mu$, we now determine their stability. The stability of a given equilibrium
is determined by small perturbations. For an equilibrium solution
$\zeta_{{\mathrm{eqm}}}(x)$, we perturb it and assume a solution of the form

$\displaystyle\zeta_{n}(x,t)=\zeta_{{\mathrm{eqm}}}(x)+\epsilon\zeta_{p}(x)e^{\lambda
t}$ (S17)
$\displaystyle\tau_{n}(t)=\tau_{{\mathrm{eqm}}}+\epsilon\tau_{p}e^{\lambda t}$
(S18)

for some unknown growth rate $\lambda$, and unknown perturbation quantities
$\zeta_{p}(x)$ and $\tau_{p}$.

Substituting this solution back into eqn. (S5)-(S7) and considering only terms
of $O(\epsilon)$, we obtain

$\displaystyle\frac{d^{4}\zeta_{p}}{dx^{4}}+\tau_{{\mathrm{eqm}}}^{2}\frac{d^{2}\zeta_{p}}{dx^{2}}+\lambda^{2}\zeta_{p}+2\tau_{p}\tau_{{\mathrm{eqm}}}\frac{d^{2}\zeta_{{\mathrm{eqm}}}}{dx^{2}}=0$
(S19)
$\displaystyle\int_{0}^{1}\frac{d\zeta_{{\mathrm{eqm}}}}{dx}\frac{d\zeta_{p}}{dx}dx=0$
(S20)
$\displaystyle\zeta_{p}(0,t)=\zeta_{p}(1,t)=\frac{d\zeta_{p}}{dx}(0,t)=\frac{d\zeta_{p}}{dx}(1,t)=0$
(S21)

The general solution of this system can be expressed as

$\displaystyle\zeta_{p}=A\cos{\mu_{+}x}+B\sin{\mu_{+}x}+C\cos{\mu_{-}x}+D\sin{\mu_{-}x}-2\tau_{p}\tau_{{\mathrm{eqm}}}\frac{\partial^{2}\zeta_{{\mathrm{eqm}}}}{\partial
x^{2}}\frac{1}{\lambda^{2}}$ (S22)
$\displaystyle\mu_{\pm}=\left(\frac{\tau_{{\mathrm{eqm}}}^{2}\pm\sqrt{{\tau_{{\mathrm{eqm}}}^{4}-4\lambda^{2}}}}{2}\right)^{1/2}.$
(S23)

Any solution must also satisfy the boundary conditions eqn. (S21).
Substituting the expressions for $\zeta_{{\mathrm{eqm}}}$ and
$\tau_{{\mathrm{eqm}}}$ into eqn. (S22) and (S21) we obtain a system of five
linear equations in five unknowns ($A,B,C,D$ and $\tau_{p}$). We can then find
(numerically) the values of $\lambda^{2}$ for which non-trivial solutions of
this homogeneous linear system exist.

Recalling from eqn. (S17) that the perturbation to the equilibrium state grows
with time like $e^{\lambda t}$, we note that unstable solutions correspond to
$\lambda^{2}>0$ while stable solutions correspond to $\lambda^{2}<0$, (i.e.
undamped oscillations). We use the sign of $\lambda^{2}$ to determine the
stability of the different equilibrium solutions already determined. The
results of the stability analysis are plotted in fig. S1c where the dashed
line indicates unstable mode and the solid line indicates that the equilibrium
is stable.

##  III Numerical simulation

In this section, we discuss the numerical simulation method used to study the
snap-through dynamics of elastic arches. Our method is based on the Discrete
Elastic Rods algorithm [2, 3].

The configuration of a planar elastic strip can be described based on its
centerline, $s$. From the numerical point of view, we discretize the strip
centerline into $N$ nodes so that the degrees of freedom are

$\mathbf{q}\equiv\left[\mathbf{x}_{0},\mathbf{x}_{1},...,\mathbf{x}_{N-1}\right]\in\mathbb{R}^{2N\times
1},$ (S24)

and the $i$-th edge vector is given by,

$\mathbf{e}_{i}=\mathbf{x}_{i+1}-\mathbf{x}_{i},\;\mathrm{with}\;i\in[0,\ldots,N-2].$
(S25)

The energy of a strip of linearly elastic material is broken down into its
stretching and bending energies. For a strip of length $L$, cross-sectional
area $w\times h$, and Young’s modulus $E$, the discrete stretching energy of
the edge connecting $\mathbf{x}_{i}$ and $\mathbf{x}_{i+1}$ is

$E^{s}_{i}=\frac{1}{2}\;EA\;\varepsilon_{i}^{2}\;||\bar{\mathbf{e}}_{i}||,$
(S26)

where $EA$ is the stretching stiffness, $||\bar{\mathbf{e}}^{i}||$ is the
undeformed length of $i$-th edge vector (quantities with a bar on top indicate
that quantity in its undeformed configuration), and $\varepsilon_{i}$ is its
axial stretch strain,

$\varepsilon_{i}=\frac{||{\mathbf{e}}_{i}||}{||\bar{\mathbf{e}}_{i}||}-1.$
(S27)

The discrete bending energy is associated with discrete curvature, and can be
written

$E^{b}_{i}=\frac{1}{2}\;EI\;\left(\kappa_{i}-\bar{\kappa}_{i}\right)^{2}\Delta\bar{L}_{i}$
(S28)

where $EI$ is the bending stiffness,
$\Delta\bar{L}_{i}={(||{\bar{\mathbf{e}}}_{i-1}||+||{\bar{\mathbf{e}}}_{i}||)}/{2}$
is the mean length of the $(i-1)$-st and $i$-th edges. Moreover, $\kappa_{i}$
is the discrete curvature,

$\kappa_{i}=\frac{2\tan({\phi_{i}}/{2})}{\Delta\bar{L}_{i}},$ (S29)

and $\phi_{i}$ is the turning angle between the two consecutive edges,
$\\{\mathbf{e}_{i-1},\mathbf{e}_{i}\\}$.

The total stretching energy of the elastic strip system can be obtained simply
by summing over all the nodes and edges, i.e.,

$E^{s}=\sum_{i=0}^{N-2}E^{s}_{i}\;\;\mathrm{and}\;\;E^{b}=\sum_{i=1}^{N-2}E^{b}_{i}.$
(S30)

It is worth noting that the strip is effectively inextensible during the snap-
through process: the elastic stretching energy penalty is so large that in
practice the uniaxial strain is zero, i.e. $\varepsilon\equiv 0$.

We numerically solve the discrete equations of motion step by step by using a
second-order symplectic Newmark-beta method [4]. At the $k$-th time step,
$t_{k}$, the generalized displacement vector $\mathbf{q}$ can be updated from
$t=t_{k}$ to $t=t_{k+1}$ by imposing the equation of motion, as

$\displaystyle\mathbb{M}\left(\Delta\mathbf{q}_{k+1}-\dot{\mathbf{q}}_{k}\;\delta
t\right)$ $\displaystyle-\frac{1}{4}\;\delta
t^{2}\;\left(\mathbf{F}^{\text{int}}_{k+1}+\mathbf{F}_{k+1}^{\text{ext}}\right)-\frac{1}{4}\;\delta
t^{2}\;\left(\mathbf{F}^{\text{int}}_{k}+\mathbf{F}_{k}^{\text{ext}}\right)=\mathbf{0}$
(S31) $\displaystyle\mathbf{q}_{k+1}$
$\displaystyle=\mathbf{q}_{k}+\Delta\mathbf{q}_{k+1}$
$\displaystyle\dot{\mathbf{q}}_{k+1}$
$\displaystyle=\frac{2}{h}\Delta\mathbf{q}_{k+1}-\dot{\mathbf{q}}_{k},$

where $\mathbb{M}$ is the diagonal mass matrix, $\delta t$ is the time step
size, $\mathbf{F}^{\text{int}}=-\nabla(E^{s}+E^{b})$ is the internal elastic
force vector, and $\mathbf{F}^{\text{ext}}$ is the external force vector
(e.g., gravity and damping). The subscript indicates the evaluation at
$t=t_{k}$ (or $t=t_{k+1}$) time step. The classical Newton-Raphson method is
applied to iteratively solve the nonlinear equations of motion. While the
error remains larger than the tolerance, we repeat the iteration process,
stopping only once the error falls below the tolerance and then moving on to
the next time step.

To achieve a clamped–clamped boundary condition, the first two nodes,
$\\{\mathbf{x}_{0},\mathbf{x}_{1}\\}$, as well as the last two nodes,
$\\{\mathbf{x}_{N-2},\mathbf{x}_{N-1}\\}$, are constrained based on the
boundary and loading conditions; all other nodes are free to evolve based on
the dynamic governing equations. The geometric and physical parameters are
identical to the experimental setup. A convergence study for both the spatial
and temporal discretizations is presented in Fig. S2, demonstrating that
$N=100$ grid points with $\delta t=0.01\mathrm{~{}ms}$ is sufficient. It is
worth noting that a small phase difference can be found in Fig. S2B. This is
caused by the implementation of the clamped boundary conditions, which means
that the effective length of the beam is slightly different for different
numbers of nodes, $N$. As a result, the natural frequencies are also slightly
different: an effect that is manifested as a small variation in the phase of
oscillations as time goes by. Nevertheless, the overall dynamic processes
(e.g., amplification) show no variations as $N$ varies.

Figure S2: Convergence study for both the temporal and spatial
discretizations. (a) Central angle $\theta$ as a function of boundary
stretching distance $-t\dot{\Delta{L}}$, with different time steps: $\delta
t\in\\{0.08,0.04,0.02,0.01\\}$ ms and a fixed number of nodes ($N=100$). (b)
Central angle $\theta$ as a function of boundary stretching distance
$-t\dot{\Delta{L}}$, with a different fixed number of nodes
($N\in\\{50,100,200,300\\}$) and a fixed time step size $\delta t=0.01$ ms.

##  IV A Toy model

Characterizing the behavior of the Elastica’s snap-through analytically is
complex because the dynamic equations involve both spatial and temporal
dependence: we would need to solve for the continuous shape of the arch at any
instant in time. To make analytical progress, we therefore choose to study a
simpler system, the double-mass von Mises truss (see Fig. 3a of the main
text). The properties of this system can be tuned so that it has the same
bifurcation behavior as the equilibrium Elastica (see Fig. 1c of the main
text), though its shape is fully described by the position of the two elements
(or masses). This model is much simpler to treat mathematically, but retains
sufficient degrees of freedom to capture both symmetric and asymmetric
transitions dynamically. This allows us to understand the key physics that
governs the amplification of the asymmetric mode.

###  IV.1 The double-mass von Mises truss Model

The double-mass von Mises truss model discussed in the main text is shown in
detail in fig. S3. It is designed in the following way: two springs of natural
length $l_{0}$ and stiffness $k$ connect two point-particles, each of mass
$m$. We refer to these as the left or right mass and spring and use a
subscript $L$ or $R$ to identify them. The two masses are connected by a third
spring of the same stiffness $k$ and natural length $l_{0C}$.

All springs are connected to the mass and to each other via torsional springs
of stiffness $B$. The torsional springs that connect the left and right
springs to the edges of the structure are relaxed when the spring makes an
angle $\alpha$ below the horizontal. This condition mimics the boundary
condition on the Elastica, in which its tangent is clamped at an angle
$\alpha$ above the horizontal. Conversely, the torsion springs connecting the
left and right springs to the central one are relaxed when the springs are
aligned or straight (i.e. when the angle between them is $\pi$) providing a
bending stiffness to the system, analogous to that of an Elastica that is
naturally flat. We therefore refer to $B$ as the bending stiffness. Without
the bending stiffness provided by the second set of torsion springs, the
system is always unstable to asymmetric deformations. Since we want to capture
the qualitative behavior of the Elastica, we need the system to be stable in
the inverted state for some values of the compression $\delta l$ and the
torsion springs are therefore an essential ingredient in the model.

Figure S3: a) Sketch of the relaxed double-mass von Mises truss with two
elements (masses) and torsion springs acting about the angles $\theta_{L}$,
$\theta_{R}$,$\theta_{LC}$ and $\theta_{RC}$. In this relaxed states, the
linear springs take their natural lengths, that is $l_{0}$ for the edge
springs and $l_{0c}$ for the edge springs. b) An example of the inverted state
of the double-mass von Mises truss after compression by $\delta l$.

##  V Static analysis of the double-mass von Mises truss

We begin by studying the stability of the inverted state of the double-mass
von Mises truss and determine its bifurcation diagram to gain insight into its
behavior. Later, we shall address the dynamics of the system and study how the
unloading rate affects the snap-through instability.

###  V.1 The energy of the double truss

In this section, we aim to write down the static (or potential) energy of the
system and study the stability of the inverted state. As shown in fig. S3, we
define the reference length of the system $l=2l_{0}+l_{0C}$ to be the sum of
the relaxed lengths of each spring. We also assume the end distance can be
changed by an amount $\delta l$, causing a deformation. In the deformed state,
we let $l_{R}$ and $l_{L}$ be the length of the left and right springs
respectively, $w_{L}$ and $w_{R}$ the horizontal distance of the left mass and
right mass from their respective edges, $\Delta h$ the difference in height
between the two masses, and $l_{C}$ and $w_{C}$ be the length of the central
spring and its horizontal width. Moreover, we let $\theta_{L}$ and let
$\theta_{R}$ be the angles the side springs make with the horizontal while
$\theta_{LC}$ and $\theta_{RC}$ are the angles between the side springs and
the central one, with $\Delta\theta$ being the angle the central spring makes
with the horizontal, as shown schematically in figure S3b). Using elementary
geometry, we have that:

$\displaystyle l_{L}$ $\displaystyle=\frac{w_{L}}{\cos\theta_{L}},$ (S32)
$\displaystyle l_{R}$ $\displaystyle=\frac{w_{R}}{\cos\theta_{R}},$ (S33)
$\displaystyle w_{C}$ $\displaystyle=l-\delta l-w_{L}-w_{R},$ (S34)
$\displaystyle\Delta h=$ $\displaystyle
w_{L}\tan\theta_{L}-w_{R}\tan\theta_{R},$ (S35) $\displaystyle l_{C}$
$\displaystyle=\sqrt{\Delta h^{2}+w_{c}^{2}},$ (S36)
$\displaystyle\Delta\theta$ $\displaystyle=\arctan(w_{c},\Delta h),$ (S37)
$\displaystyle\theta_{LC}$
$\displaystyle=(\pi/2-\theta_{L})+(\pi/2-\Delta\theta),$ (S38)
$\displaystyle\theta_{RC}$
$\displaystyle=(\pi/2-\theta_{R})+(\pi/2+\Delta\theta).$ (S39)

We can use these relations to write the energy stored in the linear springs
as:

$U_{springs}={\textstyle\frac{1}{2}}k\left[(l_{L}-l_{0})^{2}+(l_{R}-l_{0})^{2}+(l_{C}-l_{0C})^{2}\right]$
(S40)

and the bending energy stored in the torsion springs as

$U_{B}={\textstyle\frac{1}{2}}B\left[(\theta_{L}-\alpha)^{2}+(\theta_{R}-\alpha)^{2}+(\theta_{LC}-\pi)^{2}+(\theta_{RC}-\pi)^{2}\right].$
(S41)

The total energy of the system is $U=U_{springs}+U_{B}$ and can be written as
a function of the variables $w_{L}$, $w_{R}$, $\theta_{L}$, $\theta_{R}$ as
well as the control parameter $\delta l$:

$\displaystyle U=$ $\displaystyle\frac{1}{2}k\left[l_{0C}-\sqrt{\left(l-\delta
l-w_{L}-w_{R}\right)^{2}+\left(w_{L}\tan\theta_{L}-w_{R}\tan\theta_{R}\right)^{2}}\right]^{2}$
$\displaystyle+\frac{1}{2}k\left(\left(l_{0}-w_{L}\sec\theta_{L}\right)^{2}+\left(l_{0}-w_{R}\sec\theta_{R}\right){}^{2}\right)$
$\displaystyle+\frac{1}{2}B\left[\arctan\left(l-\delta
l-w_{L}-w_{R},w_{L}\tan\theta_{L}-w_{R}\tan\theta_{R}\right)+\theta_{L}\right]^{2}$
$\displaystyle+\frac{1}{2}B\left[\arctan\left(l-\delta
l-w_{L}-w_{R},w_{L}\tan\theta_{L}-w_{R}\tan\theta_{R}\right)-\theta_{R}\right]^{2}$
$\displaystyle+\frac{1}{2}B\left[\left(\alpha-\theta_{L}\right)^{2}+\left(\alpha-\theta_{R}\right)^{2}\right].$
(S42)

(Note that here $\arctan(x,y)$ is the two-argument arctangent; that is, for
$x$, $y$ in the correct quadrant, $\arctan(x,y)=\arctan(y/x)$.) To make
progress, we begin by using a suitable change of variable to emphasize the
asymmetric behavior. In particular, we introduce the symmetric and asymmetric
parts of the angle:

$\theta_{S}={\textstyle\frac{1}{2}}(\theta_{L}+\theta_{R})$

and

$\theta_{A}=\theta_{L}-\theta_{R},$

respectively. Similarly, we can decompose the change in width into symmetric,
$w_{S}=w_{L}+w_{R}-(2l_{0}-\delta l)$, and asymmetric, $w_{A}=w_{R}-w_{L}$,
parts. In this way, we may write both the angles and change in width as:

$\displaystyle\theta_{L}$
$\displaystyle=\theta_{S}+{\textstyle\frac{1}{2}}\theta_{A}$ (S43)
$\displaystyle\theta_{R}$
$\displaystyle=\theta_{S}-{\textstyle\frac{1}{2}}\theta_{A}$ (S44)
$\displaystyle w_{L}$ $\displaystyle=l_{0}-{\textstyle\frac{1}{2}}(\delta
l-w_{S}+w_{A})$ (S45) $\displaystyle w_{R}$
$\displaystyle=l_{0}-{\textstyle\frac{1}{2}}(\delta l-w_{S}-w_{A}).$ (S46)

Our primary focus in the main text is on the evolution of the angles
$\theta_{S}$ and $\theta_{A}$. To simplify this analysis and make analytical
progress, we focus now on the behavior of the energy for small values of the
clamp angle $\alpha$.

###  V.2 Energy expansion for small $\alpha$

For small $\alpha$, both the natural and inverted configurations are, in some
sense, close to the flat case (in which the masses are aligned with the edges
of the structure). This is equivalent to the small vertical displacement and
small gradient approximations used to determine the bifurcation behavior of
the continuous Elastica problem in §II and is hence also consistent with the
small values of $\alpha$ used in experiments.

We observe that for the system to be almost flat, we require the angles
$\theta_{R}$, $\theta_{L}$, as well as $\Delta\theta$, to be small. Clearly,
$\theta_{R}$, $\theta_{L}$ will be of the same size as $\alpha$. Thus, the
angle $\alpha$ determines how close we are to the flat case and is our
expansion parameter.

To understand how quantities such as $w_{R}$ and $w_{L}$ scale with $\alpha$,
let us consider the unconstrained system in which the length $l$ is free,
$\alpha$ is small, and the system can relax. In this case,
$\theta_{i}\sim\alpha$ for $i=L,R$. If we do not constrain the end-to-end
displacement, the springs will be relaxed so that $l_{i}=l_{0}$ and the width
changes according to $w_{i}\sim l_{0}\cos\theta_{i}\sim
l_{0}(1-{\textstyle\frac{1}{2}}\alpha^{2})$, meaning that the change in
horizontal length scales like $\sim\alpha^{2}$. We can thus rescale angles by
$\alpha$ and changes in horizontal lengths by $\alpha^{2}$.

We summarize our rescalings as follows:

$\displaystyle\theta_{L}$
$\displaystyle=\alpha(\psi_{S}+{\textstyle\frac{1}{2}}\psi_{A})$
$\displaystyle\theta_{R}$
$\displaystyle=\alpha(\psi_{S}-{\textstyle\frac{1}{2}}\psi_{A})$
$\displaystyle w_{L}$
$\displaystyle=l_{0}\bigl{[}1-{\textstyle\frac{1}{2}}\alpha^{2}(\mu-
W_{S}+W_{A})\bigr{]}$ $\displaystyle w_{R}$
$\displaystyle=l_{0}\bigl{[}1-{\textstyle\frac{1}{2}}\alpha^{2}(\mu-
W_{S}-W_{A})\bigr{]}$ $\displaystyle l_{0C}$ $\displaystyle=l_{0}L_{0C}$
$\displaystyle B$
$\displaystyle={\textstyle\frac{1}{3}}\beta\alpha^{2}kl_{0}^{2}$ (S47)

where $\mu=\delta l/(\alpha^{2}l)$, $\psi_{A}=\theta_{A}/\alpha$,
$\psi_{S}=\theta_{S}/\alpha$, $W_{S}=w_{S}/(l_{0}\alpha^{2})$,
$W_{A}=w_{A}/(l_{0}\alpha^{2})$. (The bifurcation parameter $\mu=\delta
l/(\alpha^{2}l)$ has emerged again, as in the continuous case.)

Note that the energy of the springs scales like the change in length squared,
meaning $U_{spring}\sim k(l_{0}\alpha^{2})^{2}$, while the bending energy
scales like the change in angle squared $U_{B}\sim B\alpha^{2}$. For the two
to balance, we need $B\sim\beta\alpha^{2}kl_{0}^{2}$ (the factor of one-third
being included to simplify later equations).

We can now expand the energy to leading order ($\alpha^{4}$). The rescaled
energy $\tilde{U}=\frac{U}{kl_{0}^{2}\alpha^{4}}$ is then given by:

$\displaystyle 192\tilde{U}=$
$\displaystyle\frac{8\psi_{A}^{2}(4\beta(L_{0C}(L_{0C}+2)+2)-3\mu
L_{0C}(5L_{0C}+4))}{L_{0C}^{2}}+\frac{3\left(L_{0C}^{2}+8\right)\psi_{A}^{4}}{L_{0C}^{2}}+48\psi_{S}^{4}$
$\displaystyle+\frac{24W_{S}\left((L_{0C}-4)\psi_{A}^{2}+4L_{0C}\left(\mu+2\mu
L_{0C}+\psi_{S}^{2}\right)\right)}{L_{0C}}+8\psi_{S}^{2}\left(9\psi_{A}^{2}+16\beta-12\mu\right)$
$\displaystyle-96W_{A}\psi_{A}\psi_{S}+48W_{A}^{2}+64\beta+48\mu^{2}(2L_{0C}(L_{0C}+2)+3)-128\beta\psi_{S}+144W_{S}^{2}$
(S48)

This energy is easily minimized with respect to variations in $W_{S}$ and
$W_{A}$ by solving $\partial U/\partial W_{A}=0$ and $\partial U/\partial
W_{S}=0$, leading to:

$\displaystyle W_{S}$
$\displaystyle=-\frac{(L_{0C}-4)\psi_{A}^{2}+4L_{0C}\left(\mu(1+2L_{0C})+\psi_{S}^{2}\right)}{12L_{0C}},$
(S49) $\displaystyle W_{A}$ $\displaystyle=\psi_{A}\psi_{S}.$ (S50)

Substituting these results back into eqn. (S48) we obtain the rescaled energy:

$\displaystyle\tilde{U}=U_{0}-a_{1}\psi_{S}+{\textstyle\frac{1}{2}}(a_{2}-a_{3}\mu)\psi_{S}^{2}+{\textstyle\frac{1}{4}}a_{4}\psi_{S}^{4}+{\textstyle\frac{1}{2}}(b_{1}-b_{2}\mu+b_{3}\psi_{S}^{2})\psi_{A}^{2}+{\textstyle\frac{1}{4}}b_{4}\psi_{A}^{4},$
(S51)

where:

$\displaystyle U_{0}$
$\displaystyle=\frac{1}{6}\left(2\beta+(2+L_{0C})^{2}\mu^{2}\right)$ (S52)
$\displaystyle a_{1}$ $\displaystyle=\frac{2}{3}\beta,\quad
a_{2}=\frac{4}{3}\beta,\quad a_{3}=2(2+L_{0C})/3,\quad a_{4}=\frac{2}{3},$
(S53) $\displaystyle b_{1}$
$\displaystyle=\frac{2+L_{0C}(2+L_{0C})}{3L_{0c}^{2}}\beta$ (S54)
$\displaystyle b_{2}$ $\displaystyle=\frac{(2+L_{0C})^{2}}{6L_{0C}},\quad
b_{3}=\frac{2+L_{0C}}{6L_{0C}},$ (S55) $\displaystyle b_{4}$
$\displaystyle=\frac{(2+L_{0C})^{2}}{24L_{0C}^{2}}.$ (S56)

Note that the signs of terms in (S51) have been chosen so that all of the
above constants are positive.

The energy of (S51) is precisely the one that appears in equation (3) of the
main text, where we have collected the prefactor of the $\psi_{A}^{2}$ term
and introduced

$f(\psi_{S};\mu):=b_{1}-b_{2}\mu+b_{3}\psi_{S}^{2}.$ (S57)

The advantage of writing the energy in this form is that its structure becomes
simple and compact, allowing us to study the equilibria of the system. From
now on, perhaps confusingly, we shall write some results with respect to the
constant defined above, and others explicitly as functions of the parameters
in the problem. This is because sometimes the results are too complex to write
explicitly as functions of the parameters, while in other cases they are
simple enough allowing us to spell them out without having to use the
constants.

###  V.3 The energy minima and their stability

The energy in equation (S51) is a simple quartic polynomial in both $\psi_{S}$
and $\psi_{A}$. To study the equilibria of the system at a given value of
$\mu$, we differentiate (S51) with respect to $\psi_{S}$ and $\psi_{A}$ and
simultaneously solve

$\frac{\partial\tilde{U}}{\partial\psi_{S}}=\frac{\partial\tilde{U}}{\partial\psi_{A}}=0.$
(S58)

####  V.3.1 Symmetric equilibria: $\psi_{A}=0$

It is immediately apparent that $\psi_{A}=0$ is one solution of the second
equilibrium equation (since $\tilde{U}$ contains only even powers of
$\psi_{A}$). To find the associated value of $\psi_{S}$ in this symmetric
state, we must solve:

$\displaystyle\frac{\partial\tilde{U}}{\partial\psi_{S}}\bigg{|}_{\psi_{A}=0}=-a_{1}+(a_{2}-a_{3}\mu)\psi_{S}+a_{4}\psi_{S}^{3}=0.$
(S59)

This equation has three zeroes (equilibria) which may not be real, depending
on the value of $\mu$. They are given by:

$\displaystyle\psi_{N}=$ $\displaystyle\frac{\sqrt[3]{2}\left(a_{3}\mu-
a_{2}\right)}{\sqrt[3]{\sqrt{108a_{4}^{3}\left(a_{2}-a_{3}\mu\right){}^{3}+729a_{1}^{2}a_{4}^{4}}+27a_{1}a_{4}^{2}}}$
$\displaystyle+\frac{\sqrt[3]{\sqrt{108a_{4}^{3}\left(a_{2}-a_{3}\mu\right){}^{3}+729a_{1}^{2}a_{4}^{4}}+27a_{1}a_{4}^{2}}}{3\sqrt[3]{2}a_{4}}$
(S60) $\displaystyle\psi_{I}=$
$\displaystyle\frac{\left(1+i\sqrt{3}\right)\left(a_{2}-a_{3}\mu\right)}{2^{2/3}\sqrt[3]{\sqrt{108a_{4}^{3}\left(a_{2}-a_{3}\mu\right){}^{3}+729a_{1}^{2}a_{4}^{4}}+27a_{1}a_{4}^{2}}}$
$\displaystyle-\frac{\left(1-i\sqrt{3}\right)\sqrt[3]{\sqrt{108a_{4}^{3}\left(a_{2}-a_{3}\mu\right){}^{3}+729a_{1}^{2}a_{4}^{4}}+27a_{1}a_{4}^{2}}}{6\sqrt[3]{2}a_{4}}$
$\displaystyle\psi_{U}$
$\displaystyle=\frac{\left(1-i\sqrt{3}\right)\left(a_{2}-a_{3}\mu\right)}{2^{2/3}\sqrt[3]{\sqrt{108a_{4}^{3}\left(a_{2}-a_{3}\mu\right){}^{3}+729a_{1}^{2}a_{4}^{4}}+27a_{1}a_{4}^{2}}}$
(S61)
$\displaystyle-\frac{\left(1+i\sqrt{3}\right)\sqrt[3]{\sqrt{108a_{4}^{3}\left(a_{2}-a_{3}\mu\right){}^{3}+729a_{1}^{2}a_{4}^{4}}+27a_{1}a_{4}^{2}}}{6\sqrt[3]{2}a_{4}}$
(S62)

These solutions correspond to the natural state $\psi_{N}$, the inverted state
$\psi_{I}$ and an unstable state in between the two, $\psi_{U}$. The inverted
and middle state annihilate when

$\frac{d\tilde{U}}{d\psi_{S}}\bigg{|}_{\psi_{A}=0}=\frac{d^{2}\tilde{U}}{d\psi_{S}^{2}}\bigg{|}_{\psi_{A}=0}=0,$

corresponding to a saddle-node bifurcation at the critical values

$\displaystyle\psi^{*}_{S}=-\left(\frac{a_{1}}{2a_{4}}\right)^{1/3}=-\left(\frac{\beta}{2}\right)^{1/3},\quad{\mu_{\ast}}^{\mathrm{vM}}=\frac{a_{2}}{a_{3}}+\frac{3}{2a_{3}}\left(2a_{1}^{2}a_{4}\right)^{1/3}=\frac{3{\psi_{S}^{*}}^{2}+2\beta}{2+L_{0C}}.$
(S63)

This means that when $\mu<{\mu_{\ast}}^{\mathrm{vM}}$, only one symmetric
equilibrium exists, the natural state $\psi_{N}$.

####  V.3.2 Asymmetric equilibria: $\psi_{A}\neq 0$

There are two other equilibrium solutions with $\psi_{A}\neq 0$. Again, since
$\tilde{U}$ involves only even powers of $\psi_{A}$, we note that if
$(\psi_{S},\psi_{A})=(a,b)$ is an equilibrium, then so is the point
$(\psi_{S},\psi_{A})=(a,-b)$. These other equilibrium points are found by
solving eqns (S58) and are best written explicitly as functions of the
parameters:

$\displaystyle\psi_{S}$
$\displaystyle=-\frac{1}{4}L_{0C}\left(L_{0C}+2\right)$ (S64)
$\displaystyle\psi_{A}$
$\displaystyle=\pm\sqrt{\frac{16L_{0C}(2+L_{0C})\mu-32(2+L_{0C}(2+L_{0C})\beta-
L_{0C}^{3}(2-L_{0C})^{3}}{2(2+L_{0C})}}$ (S65)

When the two asymmetric (saddle) equilibria meet in a pitchfork bifurcation,
locally the system ceases to be stable to asymmetric perturbations. Note that
this corresponds to the point where $f(\psi_{S};\mu)$ changes from being
positive to being negative. This occurs at a critical value,
$\mu=\mu_{A}^{\mathrm{vM}}$, that can be calculated by solving
${\textstyle\frac{1}{2}}(b_{1}+b_{2}\mu+b_{3}\psi_{I}^{2})=0$, telling us that
the pitchfork bifurcation occurs at
$\psi_{S}=-\frac{1}{4}L_{0C}\left(L_{0C}+2\right)$ and at the critical value
of

$\displaystyle\mu_{A}^{\mathrm{vM}}=\frac{1}{16}L_{0C}^{2}\left(2+L_{0C}\right)^{2}+\frac{2(2+L_{0C}(2+L_{0C})\beta}{L_{0C}(2+L_{0C})^{2}}.$
(S66)

####  V.3.3 Reconstructing the elastica’s bifurcation diagram

Our static analysis has revealed that the structure of the energy, and thus
the behavior of the double-mass von Mises truss, is entirely determined by the
choice of three dimensionless parameters. These parameters are $\beta$,
indicating the ratio between bending and stretching stiffness, $L_{0C}$,
representing the length ratio between the side and central springs, and $\mu$,
which characterizes the compression in the system. The latter parameter can be
easily controlled and is used to induce the snap-through in the system.
Conversely, the first two are material parameters and are kept fixed.

To run all our numerical calculations in the paper, we choose $\beta$ and
$L_{0C}$ so that the bifurcation diagram of the double-mass von Mises truss is
as close as possible to that of the shallow arch derived in §II. In
particular, we require that the saddle-node bifurcation point described in
equation (S63) occurs at $\mu_{*}=0.247$ (i.e. we ensure that
${\mu_{\ast}}^{\mathrm{vM}}={\mu_{\ast}}$) and that the loss of stability to
asymmetric modes found in equation $\eqref{eqn:pitchfork}$ occurs at
$\mu_{A}^{\mathrm{vM}}=\mu_{A}=1/4$. These two conditions lead to a choice of
$\beta=0.110$ and $L_{0C}=0.654$.

###  V.4 Stability of the inverted state

If we consider the system in the inverted state with $\psi_{A}=0$ and imagine
slowly decreasing $\mu$ from some large value, we may then ask the question:
_when and how does the system go unstable?_ To answer, we begin by only
considering the quasi-static problem and thus ignore dynamic effects.

At first, the system is in equilibrium, with $\psi_{S}=\psi_{I}$. However, at
$\mu=\mu_{A}$, the coefficient of $\psi_{A}^{2}$ in the energy becomes
negative. This occurs once the inverted state enters the amplification region
(discussed in the main text) given by

$f(\mu,\psi_{S})<0\,\,\Rightarrow\,\,\,\psi_{S}^{2}<\frac{b_{2}\mu-
b_{1}}{b_{3}},$ (S67)

and thus becomes unstable to asymmetric modes. At this point, the system snaps
through via asymmetric modes to the only stable equilibrium (minimum): the
natural state. Therefore, a quasi-static analysis suggests that our double-
mass von Mises truss becomes unstable to asymmetric modes first, like an
Elastica.

If we consider dynamic effects, however, the behavior changes. If the end
distance is changed fast, then the (dynamic) value of $\psi_{S}$ lags behind
the actual value of the equilibrium $\psi_{I}$ due to inertial effects. This
means that $\psi_{S}^{2}>\psi_{I}^{2}$, resulting in a larger value of the
coefficient $f(\psi;\mu)$ shown in eqn (S57). If the dynamics is sufficiently
fast, the system lags behind enough to never enter the amplification region
and thus never becomes unstable to asymmetric modes. However, a (delayed)
snap-through still occurs when $\mu<{\mu_{\ast}}$, as the inverted state
ceases to exist altogether in a saddle-node bifurcation and the system is
forced to jump to the natural state.

Therefore, dynamic effects determine whether the system enters or not the
amplification region and thus characterise both the symmetry and time of the
snap-through. To make this claim more precise, we need to derive the equations
of motion and solve for the path of the system.

##  VI Dynamics of the double-mass von Mises truss

In this section, we supplement the elastic energy already considered with the
kinetic energy of the system. This allows us to study the dynamics and thereby
elucidate how the rate $|\dot{\mu}|$ at which the end-to-end displacement is
decreased modifies the behavior of the system. Here, we show that for
sufficiently small $|\dot{\mu}|$ (almost quasi-static variations), the double-
mass von Mises truss snaps from the inverted to the natural state via an
asymmetric mode. However, for larger $|\dot{\mu}|$, the importance of
asymmetric modes to promote the snap through decreases until, for
$|\dot{\mu}|>|\dot{\mu}|_{c}$, the double-mass von Mises truss always goes
unstable symmetrically.

###  VI.1 Deriving the equations of motion

Since we have already derived the potential energy of the system, it is
convenient to write down the kinetic energy and then use Lagrangian mechanics
to derive the equations of motion. The kinetic energy is easily written down
by considering the motion (in polar coordinates) of the masses about the two
edge elements:

$T={\textstyle\frac{1}{2}}m\left[(l_{L}\dot{\theta}_{L})^{2}+\dot{l}_{L}^{2}+(l_{R}\dot{\theta}_{R})^{2}+\dot{l}_{R}^{2}\right]$
(S68)

where dots indicate derivatives with respect to time $t$. We can now use the
rescalings from equation (S47), while also rescaling time by the typical time
scale of elastic oscillations caused by the torsion springs (which scales like
$\sim\sqrt{B/ml_{0}^{2}}\sim\sqrt{k/m}\alpha$), i.e. we introduce the
dimensionless time $\tau$ defined by

$t=\left(\sqrt{\frac{m}{k}}\frac{1}{\alpha}\right)\tau.$ (S69)

Expanding in small $\alpha$ we can write the rescaled kinetic energy
$\tilde{T}=\frac{T}{kl_{0}^{2}\alpha^{4}}$ as

$\tilde{T}=\dot{\psi}_{S}(\tau)^{2}+{\textstyle\frac{1}{4}}\dot{\psi}_{A}^{2}+H.O.T.$
(S70)

where now, in a slight abuse of notation, dots refer to derivative with
respect to the rescaled time $\tau$. Note that the kinetic energy scales like
the potential energy in $\alpha$ and the two are thus comparable. The terms in
the kinetic energy can be explained by realizing that, for small angles, only
the vertical velocities of the masses matter. This is because the change in
vertical displacements scale like $\alpha^{2}$ (from the scaling in equation
(S47)) while vertical displacements scale like $\alpha$ giving
$y_{L}=l_{0}\sin(\alpha(\psi_{S}+{\textstyle\frac{1}{2}}\psi_{A}))\sim
l_{0}\alpha(\psi_{S}+{\textstyle\frac{1}{2}}\psi_{A})$ and
$y_{R}=l_{0}\sin(\alpha(\psi_{S}-{\textstyle\frac{1}{2}}\psi_{A}))\sim
l_{0}\alpha(\psi_{S}-{\textstyle\frac{1}{2}}\psi_{A})$. Differentiating these
positions with respect to time leads to
$\dot{y_{L}}\sim\alpha^{2}\sqrt{k/m}l_{0}(\dot{\psi_{S}}+{\textstyle\frac{1}{2}}\dot{\psi_{A}})$
and
$\dot{y_{R}}\sim\alpha^{2}\sqrt{k/m}l_{0}(\dot{\psi_{S}}-{\textstyle\frac{1}{2}}\dot{\psi_{A}})$.
Substituting in the kinetic energy
${\textstyle\frac{1}{2}}m(\dot{y}_{L}^{2}+\dot{y}_{R}^{2})$ and rescaling
yields the result in equation (S70).

We thus write the (rescaled) Lagrangian of the system $L=\tilde{T}-\tilde{U}$,
where $\tilde{U}$ is given in eqn. (S48). We minimize the action with respect
to variations in $W_{S}$, $W_{A}$, $\psi_{S}$ and $\psi_{A}$. Minimization
with respect to $W_{S}$ and $W_{A}$ leads to (S49) and (S50), and thus yields
no new information. Minimization with respect to $\psi_{A}$ and $\psi_{S}$
leads to the equations of motion:

$\displaystyle\ddot{\psi}_{S}$
$\displaystyle={\textstyle\frac{1}{2}}a_{1}-{\textstyle\frac{1}{2}}(a_{2}-a_{3}\mu)\psi_{S}-{\textstyle\frac{1}{2}}a_{4}\psi_{S}^{3}-{\textstyle\frac{1}{2}}b_{3}\psi_{S}\psi_{A}^{2}$
(S71) $\displaystyle\ddot{\psi}_{A}$
$\displaystyle=-2f(\psi_{S},\mu)\psi_{A}-2b_{4}\psi_{A}^{3}$ (S72)

where the right-hand side of both equations is just the derivative of
$\tilde{U}$ in eqn. (S51) with respect to $\psi_{S}$ and $\psi_{A}$ for the
first and second equations, respectively. For small $\psi_{A}$, these
equations coincide with eqns. (5) and (6) of the main text.

Unfortunately, these equations are non-linear and cannot be solved explicitly.
Furthermore, the system does not conserve energy: substituting an explicit
time dependence of the end-shortening on time, shows that the Lagrangian $L$
depend explicitly on time. Note that in the main text, we substitute
$\mu=\mu_{A}-|\dot{\mu}|\tau$ so that the quasi-static instability occurs at
$\tau=0$. For the numerical calculations, we change this convention and we
write $\mu=\mu_{i}-|\dot{\mu}|\tau$ where $\mu_{i}$ is our initial end-
shortening. This choice is made for convenience and should not cause concern
since all numerical results are plotted with respect to $\mu$, not $\tau$,
making any shift in time irrelevant.

We can plot the trajectories of the system numerically. However, to do so,
some care has to be taken when imposing the initial conditions.

####  VI.1.1 Initial conditions

In principle, we could choose any reasonable initial condition, as for
instance $\psi_{S}(\tau_{i})=\psi_{i}$, $\dot{\psi}_{S}(\tau_{i})=0$ and
$\psi_{A}(\tau_{i})=\psi^{0}_{A}$, $\dot{\psi}_{A}(\tau_{i})=0$. However, our
choice of initial condition will affect the trajectory of the system, thus
adding extra variables to the problem. For example, choosing an initial
$\psi_{S}$ that is far form the equlibrium inverted solution, would introduce
large oscillations that may heavily affect the dynamics. Similarly, starting
on the equilibrium branch, but with a ‘wrong’ initial $\dot{\psi}_{S}$, may
also introduce large oscillations and thus the same issue.

To get around this problem, we require the solution to the equations of motion
to map back to the stable inverted solution as $\tau\to-\infty$, i.e. at a
large time prior to the instability. (This is analogous to the method used to
find the correct initial conditions when the control parameter $\mu$ varies
dynamically to control snap-through [5]). Imposing such a condition on both
equations of motion is difficult. We can assume that the amplitude of the
asymmetric component of motion is much smaller than that of the symmetric one,
and ignore the contribution of $\psi_{A}$ for large negative times. We then
need to find an expansion for $\tau\to-\infty$ of

$\ddot{\psi}_{S}={\textstyle\frac{1}{2}}a_{1}-{\textstyle\frac{1}{2}}(a_{2}-a_{3}\mu)\psi_{S}-{\textstyle\frac{1}{2}}a_{4}\psi_{S}^{3}.$
(S73)

We expand $\psi_{S}$ in powers of $\sqrt{-\tau}$ and solve order by order.
Substituting in the definitions of the constants $a_{1}$, $a_{2}$, etc. we
find that:

$\displaystyle\psi_{S}=-\sqrt{(2+L_{0C})}\sqrt{|\dot{\mu}|}\sqrt{-\tau}+\frac{\beta}{\sqrt{(2+L_{0C})}\sqrt{|\dot{\mu}|}\sqrt{-\tau}}-\frac{\beta}{2(2+L_{0C})|\dot{\mu}|\tau}+....$
(S74)

Correspondingly, the first derivative of $\psi_{S}$ for large negative times
is:

$\displaystyle\dot{\psi}_{S}=\frac{\sqrt{(2+L_{0C})}\sqrt{|\dot{\mu}|}}{2\sqrt{-\tau}}-\frac{\beta}{2\sqrt{(2+L_{0C})}\sqrt{|\dot{\mu}|}(-\tau)^{3/2}}+\frac{\beta}{2(2+L_{0C})|\dot{\mu}|\tau^{2}}+....$
(S75)

In the numerical solutions presented in figure 4 of the main text, we use the
first two terms in the above series to determine $\psi_{S}(\tau_{i})$ and
$\dot{\psi}_{S}(\tau_{i})$ evaluated at initial time $\tau_{i}=-100$ to impose
boundary conditions on $\psi_{S}$. This choice ensures that the initial error
on $\psi_{S}$ is of order $10^{-4}$ while the error on its derivative is of
order $10^{-6}$, as evaluated from the first discarded term in the above
series. We find this approximation sufficient so that the numerical solution
correctly traces back to the equilibrium branch at large negative times
without any unwanted oscillations.

The boundary conditions for $\psi_{A}$ are considerably simpler: we just
impose $\dot{\psi}_{A}(\tau_{i})=0$ and an initial perturbation
$\psi_{A}(\tau_{i})=\psi^{0}_{A}$.

##  VII Critical path

As discussed in the main text, the initial asymmetry in the system is
amplified only when the trajectory of the system enters the amplification
region, which happens when $|\dot{\mu}|>|\dot{\mu}_{c}|$. In this section, we
determine the critical rate $|\dot{\mu}_{c}|$ that marks this transition
between trajectories that do and do not enter the amplification region.

###  VII.1 Approximation near the saddle-node bifurcation point

When the perturbation in the asymmetric mode is small, the dynamic system can
be simplified to:

$\displaystyle\ddot{\psi}_{S}$
$\displaystyle={\textstyle\frac{1}{2}}a_{1}-{\textstyle\frac{1}{2}}(a_{2}-a_{3}\mu)\psi_{S}-{\textstyle\frac{1}{2}}a_{4}\psi_{S}^{3}$
(S76) $\displaystyle\ddot{\psi}_{A}$
$\displaystyle=-2(b_{1}-b_{2}\mu+b_{3}\psi_{S}^{2})\psi_{A}$ (S77)

It is instructive to first write eqn. (S76) in terms of the underlying
parameters, $\mu$ and $\beta$; this leads to:

$3\ddot{\psi}_{S}=\beta-\left((2+L_{0C})\mu-2\beta\right)\psi_{S}-\psi_{S}^{3}.$
(S78)

We can approximate the trajectory near the saddle-node bifurcation point given
by $\psi_{S}^{*}$ and $\mu_{*}$ defined in equation (S63). To do so, we let:

$\displaystyle\psi(\mu)$
$\displaystyle=\psi^{*}_{S}+|\dot{\mu}|^{2/5}\left(\frac{2}{\beta}\right)^{1/15}\left(\frac{3}{2-L_{0C}}\right)^{2/15}\delta\psi_{S}(\Delta)$
(S79) $\displaystyle\mu$
$\displaystyle={\mu_{\ast}}+|\dot{\mu}|^{4/5}\left(\frac{2}{\beta}\right)^{2/15}\left(\frac{3}{2-L_{0C}}\right)^{1/15}\Delta$
(S80)

and keep only the leading order terms in $\delta\psi_{S}$ and $\Delta$ so that
eqn. (S76) takes the simple form:

$\displaystyle\frac{d^{2}}{d\Delta^{2}}\delta\psi_{S}$
$\displaystyle=\delta\psi_{S}^{2}(\Delta)-\Delta.$ (S81)

This equation can easily be solved numerically and describes the trajectory of
our system near the saddle-node bifurcation for all parameters $\beta$,
$L_{0C}$ and $|\dot{\mu}|$.

Next, we focus on understanding when the trajectory of the system enters the
amplification region. To do so, we look at the function $f$ in eqn. (S57). We
are interested in understanding which values of $|\dot{\mu}|$ result in $f$
changing sign and becoming negative, thus amplifying the asymmetric
perturbation. We substitute the new variables from eqn. (S80) in the
definition of $f$ and obtain

$\displaystyle
f(\Delta)=-c_{0}-c_{1}\left(3|\dot{\mu}|^{4/5}\Delta-\left(|\dot{\mu}|^{2/5}\delta\psi_{S}-c_{2}\right)^{2}\right),$
(S82)

where

$\displaystyle
c_{0}=\frac{6L_{0C}(2\beta)^{1/3}+3L_{0C}^{2}(2\beta^{2})^{1/3}-8\beta}{12L_{0C}},\quad\quad$
(S83) $\displaystyle
c_{1}=\left(\frac{(2+L_{0C})}{3}\right)^{9/5}\frac{1}{2^{13/15}L_{0C}\beta^{2/15}},\quad
c_{2}=\left(\frac{3\beta}{2(2+L-{0C})}\right)^{2}.$ (S84)

The critical path, where the trajectory only touches but never enters the
amplification region, can be found by solving $f(\Delta)=0$ and
$f^{\prime}(\Delta)=0$ simultaneously. We can easily do this numerically and
we obtain that the critical rate is $|\dot{\mu}_{c}|=0.01014$, as discussed in
the main text. When $|\dot{\mu}|<|\dot{\mu}_{c}|$, the trajectory enters the
bifurcation region and amplification of the asymmetric mode is observed.

## References

  * [1] Lowess smoothing, https://www.mathworks.com/help/curvefit/lowess-smoothing.html.
  * Bergou  _et al._ [2008] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, and E. Grinspun, Discrete elastic rods, in  _ACM Trans. Graph._ , Vol. 27 (ACM, 2008) p. 63.
  * Bergou  _et al._ [2010] M. Bergou, B. Audoly, E. Vouga, M. Wardetzky, and E. Grinspun, Discrete viscous threads, in  _ACM Trans. Graph._ , Vol. 29 (ACM, 2010) p. 116.
  * Huang and Jawed [2019] W. Huang and M. K. Jawed, Newmark-beta method in discrete elastic rods algorithm to avoid energy dissipation, J. Appl. Mech. 86, 084501 (2019).
  * Liu  _et al._ [2021] M. Liu, M. Gomez, and D. Vella, Delayed bifurcation in elastic snap-through instabilities, J. Mech. Phys. Solids 151, 104386 (2021).