# Weakly nonlinear analysis of Turing pattern dynamics on curved surfaces

Ryosuke Nishide1<EMAIL_ADDRESS>Shuji Ishihara1,2 1
Graduate School of Arts and Sciences, The University of Tokyo, Komaba 3-8-1,
Meguro-ku, Tokyo 153-8902, Japan
2 Universal Biology Institute, The University of Tokyo, Komaba 3-8-1, Meguro-
ku, Tokyo 153-8902, Japan

###### Abstract

Pattern dynamics on curved surfaces are ubiquitous. Although the effect of
surface topography on pattern dynamics has gained much interest, there is a
limited understanding of the roles of surface geometry and topology in pattern
dynamics. Recently, we reported that a static pattern on a flat plane can
become a propagating pattern on a curved surface [Nishide and Ishihara, Phys.
Rev. Lett. 2022]. By examining reaction-diffusion equations on axisymmetric
surfaces, certain conditions for the onset of pattern propagation were
determined. However, this analysis was limited by the assumption that the
pattern propagates at a constant speed. Here, we investigate the pattern
propagation driven by surface curvature using weakly nonlinear analysis, which
enables a more comprehensive approach to the aforementioned problem. The
analysis reveals consistent conditions of the pattern propagation similar to
our previous results, and further predicts that rich dynamics other than
pattern propagation, such as periodic and chaotic behaviors, can arise
depending on the surface geometry. This study provides a new perspective on
the relationship between surfaces and pattern dynamics and a basis for
controlling pattern dynamics on surfaces.

##  I Introduction

Pattern formation occurs at a wide range of scales in natural and engineered
systems, including nanoscale metal stripes, chemical waves, fish and animal
skin, vegetation, and atmospheric convection [4, 1, 2, 5, 3, 6, 7, 8]. The
stages in which pattern formation occurs are diverse in terms of spatial
dimensions and geometries, among which curved surfaces are common and
fundamental. In particular, biological systems are rich sources of pattern
dynamics on curved surfaces [9, 10, 12, 13, 14, 15, 11]. Recent developments
in experimental techniques for measuring and controlling curved surfaces have
revealed the functional roles of curved surfaces in pattern formation [16, 17,
18]. Bin/Amphiphysin/Rvs domain proteins sense the curvature of cellular
membranes and regulate the cellular shape [19, 20]. The distribution of the
mechanosensitive protein Piezo1 is regulated by membrane curvature [21].
Cellular migration is guided by the curvature of a substrate and is called
“curvotaxis” [22]. The thickness of a cellular sheet depends on the curvature
of the substrate [23].

Several theoretical studies have revealed the effects of geometry and topology
on pattern formation and dynamics. Geodesic curvature causes the splitting
[24], rectification [25], and swirling [26, 27, 14] of excitable waves and
variations in the wavefront speed [30, 28, 29]. It also controls the position
of the localization pattern [31, 32]. The rectification of propagating
patterns by curved surfaces was also reported in the collective motion of
self-propelled particles [33]. For pattern dynamics described by polar and
nematic variables, topological defects are constrained by Poincaré–Hopf
theorem on closed surfaces [34], and have been investigated in liquid crystals
[35, 36], flocking [37], and active nematics [38]. This restriction plays a
critical role in hydra regeneration, where singular points of cortical actin
fibers dictate the positional information of morphogenesis [15].

Turing patterns, a prominent example of pattern formation [39, 40, 41, 42],
have also attracted much interest in their dynamics on curved surfaces; in
fact, Turing himself studied the pattern on a sphere to consider embryogenesis
[39]. Since then, Turing patterns have been studied on several curved surfaces
[43, 44, 46, 45]: spheres [39, 47, 48, 49, 50, 51], hemispheres [52, 53], tori
[54, 51], ellipsoids [54, 55], and deformed axisymmetric cylinders and spheres
[56]. These studies have reported that the profile and position of a Turing
pattern change to reflect the surface shape. For example, Frank et al. studied
pattern formation on a cylinder with a ridge, and found that the stripe
pattern position is modulated by the ridge, termed as “pinning” [56]. All
these studies presumed that a Turing pattern, which is static on a flat plane,
remains static irrespective of the surface geometry.

We recently reported that a Turing pattern can no longer remain static on
curved surfaces [57], indicating a new mechanism of pattern propagation caused
by surface curvature, i.e., a pattern that is static on a flat plane
propagates on a curved surface. We primarily studied patterns on axisymmetric
surfaces, which enabled a tractable analysis. Linear stability analysis
indicates no sign of Hopf bifurcation, similar to the standard Turing
instability condition, and the growth rate near the uniform state is
indistinguishable between flat planes and curved surfaces (see Sec. II.2).
Thus, analysis incorporating nonlinearity is required to understand the onset
of such time-dependent dynamics. We performed numerical simulations and
theoretical analyses based on the relationship between the propagation
velocity and the pattern profile (see Eq. (7) in Sec. II.3) and revealed that
the onset of propagation is conditioned by symmetries of the surface and
pattern. The helical pattern propagates unless the surface is reflection-
symmetric, whereas the parallel pattern is always static (see Fig. 1). These
results are generic and were confirmed by numerical simulations of the
Brusselator [58, 59] and Lengyel–Epstein [60, 61] models.

In this paper, we adopt an alternative analysis method to uncover the manner
in which a static pattern on a plane can propagate on a curved surface. Our
previous analysis assumed that patterns propagate at a constant velocity along
the azimuthal direction and used the relationship in Eq. (7). Instead, in this
study, we consider the solution of the reaction-diffusion equation (RDE) in
the vicinity of the Turing bifurcation point and perform a weakly nonlinear
analysis. Near the bifurcation point, the amplitudes of the pattern profile
are small, which allows us to expand the pattern by the superposition of a few
linear modes and derive the amplitude equations [62, 59, 63, 64, 46]. The
obtained amplitude equations indicate that different mode-mode interactions
occur depending on the symmetries of the surface and pattern profile,
resulting in either the onset or suppression of pattern propagation. This
analysis not only leads to consistent results as in [57], but further predicts
that non-trivial dynamics other than propagating patterns, including periodic
and chaotic behaviors, can appear on curved surfaces. These results are
confirmed by numerical simulations.

This paper is organized as follows. Sec. II describes the setup of our
problem, and summarizes our previous results that the Turing pattern can
propagate on a curved surface. In Sec. III, the amplitude equations for the
RDE on an axisymmetric surface are discussed. Possible forms of amplitude
equations depending on the surface and pattern symmetries are determined. In
Sec. IV, the correspondence between the steady-state solutions of the
amplitude equations and static and propagating patterns is discussed. In Sec.
V, the amplitude equations for the Brusselator model are derived and the
discussions in the previous sections are numerically confirmed. Sec. VI shows
that surface curvature can cause diverse pattern dynamics other than
propagation, which is a summary of the joint paper [65]. Finally, the results
are summarized in Sec. VII.

##  II Model and Propagating Pattern

###  II.1 Model

We study reaction-diffusion systems on an axisymmetric surface represented by
$\bm{r}=(x,r(x)\cos\theta,r(x)\sin\theta)$ (Fig. 2(a)). The coordinates $x$
and $\theta$ are defined as $-2\pi\leq x<2\pi$ and $0\leq\theta<2\pi$,
respectively, and periodic boundary conditions are used for both axes unless
otherwise mentioned. We consider the general form of the RDE as follows:

$\displaystyle\partial_{t}\bm{u}=D\Delta\bm{u}+\bm{R}(\bm{u}),$ (1)

where ${\bm{u}}$ is the vector representation of the chemical concentrations
and a function of the position on the surface $(x,\theta)$ and the time $t$,
$D$ is a diagonal matrix composed of the diffusion coefficients, $\Delta$ is
the Laplace–Beltrami operator, and $\bm{R}(\bm{u})$ is the reaction term. The
Laplace–Beltrami operator $\Delta$ characterizes the effect of the surface
curvature on the system [67]. In general, for a curved surface
$\bm{r}=\bm{r}(x^{1},x^{2})$ parameterized by the coordinates $(x^{1},x^{2})$,
the Laplace–Beltrami operator $\Delta$ for a scalar-field variable is
determined by the metric tensor
$g_{ij}=\partial_{i}\bm{r}\cdot\partial_{j}\bm{r}$, where we use the notation
$\partial_{i}\equiv\partial/\partial{x^{i}}$:

$\displaystyle\Delta\bullet=\frac{1}{\sqrt{g}}\partial_{i}\bigg{(}\sqrt{g}g^{ij}\partial_{j}\bullet\bigg{)}~{},$
(2)

where, $g^{ij}$ is the inverse of $g_{ij}$, $g\equiv|\det(g_{ij})|$, and we
use the Einstein summation convention for repeated indices. For an
axisymmetric surface, $g_{xx}=1+r^{\prime 2}$, $g_{\theta\theta}=r^{2}$,
$g_{x\theta}=g_{\theta x}=0$, and $\sqrt{g}=r\sqrt{1+r^{\prime 2}}$, where
$r^{\prime}\equiv dr(x)/dx$. The Laplace–Beltrami operator is given by:

$\displaystyle\Delta\bullet=\frac{1}{r\sqrt{1+r^{\prime
2}}}\partial_{x}\bigg{(}\frac{r}{\sqrt{1+r^{\prime
2}}}\partial_{x}\bullet\bigg{)}+\frac{1}{r^{2}}\partial^{2}_{\theta}\bullet~{}.$
(3)

We assume that Eq. (1) has a unique uniform steady solution, $\bm{u}_{0}$,
that satisfies $\bm{R}(\bm{u}_{0})=0$, and the steady-state solution becomes
unstable because of the Turing instability indicated by the dispersion
relation discussed below.

As a representative example, we use the Brusselator model given by

$\displaystyle\partial_{t}u$ $\displaystyle=D_{u}\Delta u+u^{2}v-bu-u+a~{},$
$\displaystyle\partial_{t}v$ $\displaystyle=D_{v}\Delta v-u^{2}v+bu~{},$ (4)

where $u(t,x,\theta)$ and $v(t,x,\theta)$ represent the chemical
concentrations at the position on the surface $(x,\theta)$ and time $t$,
$D_{u}$ and $D_{v}$ are the diffusion coefficients of $u$ and $v$,
respectively, and $a$ and $b$ are model parameters. The Brusselator model has
a uniform steady solution, $\bm{u}_{0}=(a,b/a)^{T}$, and, although it can show
various spatiotemporal patterns including static stripes, spiral propagation,
and oscillatory patterns [66], in this study we selected parameter sets such
that the system exhibits a Turing pattern on a flat surface. For the numerical
data presented in this study, we verified that the pattern is static on a flat
plane. The details of the numerical simulation method are provided in
Supplemental Text Sec. S1 [68].

###  II.2 Linear stability analysis

Turing instability is characterized by the dispersion relation $\mu(\lambda)$,
where $\lambda$ is the eigenvalue of the Laplace–Beltrami operator. The
dispersion relation $\mu(\lambda)$ is obtained by a linear stability analysis
at the uniform state ${\bm{u}}={\bm{u}}_{0}$, as an eigenvalue of the linear
operator

$\displaystyle\mathcal{L}=D\Delta+\partial_{\bm{u}}{\bm{R}}(\bm{u}_{0})~{},$
(5)

where $\partial_{\bm{u}}{\bm{R}}(\bm{u}_{0})$ denotes the Jacobian matrix of
${\bm{R}}(\bm{u})$ at ${\bm{u}}={\bm{u}_{0}}$. The eigenfunctions of
$\mathcal{L}$ are represented by
$\bm{U}_{\lambda,k}(x,\theta)={\bm{A}}_{\lambda}W_{\lambda,k}(x,\theta)$, with
the Laplace–Beltrami eigenfunction $W_{\lambda,k}(x,\theta)$ defined by
$\Delta W_{\lambda,k}=-\lambda W_{\lambda,k}$. For flat surfaces, $\lambda$
coincides with the square of the wavenumber. For an axisymmetric surface,
$W_{\lambda,k}$ can be factored into
$W_{\lambda,k}(x,\theta)=X_{\lambda,k}(x)e^{ik\theta}$, where $k$ is the
wavenumber along the $\theta$-direction and is an integer. $X_{\lambda,k}(x)$
is a function of $x\in[-2\pi,~{}2\pi)$ and is the solution to the following
equation:

$\displaystyle\frac{r}{\sqrt{1+r^{\prime
2}}}\partial_{x}\bigg{(}\frac{r}{\sqrt{1+r^{\prime
2}}}\partial_{x}X_{\lambda,k}\bigg{)}+\left(\lambda
r^{2}-k^{2}\right)X_{\lambda,k}=0~{},$ (6)

under the periodic boundary condition. Note that the eigenvalue $\lambda$ is a
nonnegative real number ($\lambda\in\mathbb{R}_{0}^{+}$), which is a general
feature of the Laplace–Beltrami operator. In addition, the eigenvalues are
discretized for a closed surface with a finite area. The eigenvalues are
degenerate two-fold for $k\neq 0$, for which $W_{\lambda,k}(x,\theta)$ and
$W_{\lambda,-k}(x,\theta)$ are pairwise eigenfunctions that correspond to the
anti-phase shift along the $\theta$-axis. Vector ${\bm{A}}_{\lambda}$ is the
eigenvector of the matrix $-\lambda D+\partial_{\bm{u}}{\bm{R}}(\bm{u}_{0})$,
where $\mathcal{L}{\bm{U}}_{\lambda,k}=\mu(\lambda){\bm{U}}_{\lambda,k}$ is
satisfied. The dispersion relation $\mu(\lambda)$ represents the growth rate
of a mode characterized by $\lambda$. For Turing instability, $\mu(\lambda)$
takes real positive values within a finite range of $\lambda$, and the
corresponding eigenvector $\bm{A}_{\lambda}$ is a real-valued vector. Note
that the dispersion relation $\mu(\lambda)$, and thus the Turing condition, is
the same for any surface.

###  II.3 Pattern propagation on axisymmetric surface

By setting the parameters in the Brusselator model to
$(a,b,D_{u},D_{v})=(2.0,4.5,0.5,1.8)$, we obtain the dispersion relation shown
in Fig. 2(b). In the numerical simulation of the RDE, the initial small
fluctuation around the uniform solution $\bm{u}_{0}$ grows and the chemical
concentration eventually forms a spatially periodic pattern (i.e., the Turing
pattern). The pattern is static on a plane (Fig. 2(c)), whereas it propagates
on an axisymmetric surface, as shown in Fig. 2(d). Futhermore, the propagation
is observed when the Neumann boundary condition is employed (Fig. S1 [68]).

The velocity of the pattern is related to the pattern profile and surface
shape. Assuming a traveling wave solution $\bm{u}(x,\rho)$ with
$\rho=\theta-\omega t$, the angular velocity $\omega$ along the $\theta$-axis
is obtained by the following equation:

$\displaystyle\omega=-\frac{\langle\partial_{\rho}\bm{u},\bm{R}\rangle}{\langle\partial_{\rho}\bm{u},\partial_{\rho}\bm{u}\rangle}~{},$
(7)

where the inner product is defined by

$\displaystyle\langle\bm{\Phi}(x^{1}\\!,x^{2}),\\!\bm{\Psi}(x^{1}\\!,x^{2})\rangle\equiv\\!\int\\!dS~{}\bm{\Phi}^{\dagger}(x^{1}\\!,x^{2})\bm{\Psi}(x^{1}\\!,x^{2})~{},$
(8)

with the area element $dS=\sqrt{g}dx^{1}dx^{2}$, and $\bm{\Phi}^{\dagger}$
denoting the Hermitian adjoint vector of $\bm{\Phi}$.

In our previous study [57], we used Eq. (7) to reveal the relationship between
pattern propagation and the symmetry of patterns and surfaces, as summarized
in Fig. 1. The analysis therefore largely depended on the assumption that the
pattern propagates at constant speed. Note also that we primarily studied
axisymmetric surfaces with $2\pi$-periodicity along the $x$-axis; in this
study, we will show that the loss of the $2\pi$-periodicity results in the
propagation of the parallel pattern.

##  III Amplitude equation and symmetry of system

###  III.1 Mode decomposition of the pattern dynamics

We analyze the propagating pattern on an axisymmetric surface by expanding the
pattern using eigenfunctions of the linear operator in a uniform state,
${\bm{U}}_{\lambda,k}(x,\theta)={\bm{A}}_{\lambda}W_{\lambda,k}(x,\theta)$. If
the pattern is not significantly far from the bifurcation point of the Turing
instability, where the uniform solution ${\bm{u}}_{0}$ becomes unstable, the
pattern profile can be approximated by a linear combination of a few
eigenfunctions ${\bm{U}}_{\lambda,k}$. We evaluate the expansion coefficients
by numerically calculating the inner product between the eigenfunction
${\bm{U}}_{\lambda,k}(x,\theta)$ and the propagating pattern
${\bm{U}}(t,x,\theta)\equiv{\bm{u}}(t,x,\theta)-{\bm{u}}_{0}$.

This expansion was applied to the numerically obtained Brusselator patterns
(see Figs. 3 and S2–S4 [68]). The helical propagating pattern shown in Figs.
3(a) and (b) is well approximated by two pairwise eigenfunctions with the same
wavenumber, $k=2$. Fig. 3(c) shows the main part of the expansion:

$\displaystyle\bm{U}_{\ast}(t,x,\theta)$
$\displaystyle=C_{1}(t)\bm{U}_{\lambda_{1},k}(x,\theta)+C_{2}(t)\bm{U}_{\lambda_{2},k}(x,\theta)$
(9)
$\displaystyle~{}~{}~{}+\bar{C}_{1}(t)\bm{U}_{\lambda_{1},-k}(x,\theta)+\bar{C}_{2}(t)\bm{U}_{\lambda_{2},-k}(x,\theta)~{}.$

This closely approximates the original patterns $U=u(t,x,\theta)-a$ and
$V=v(t,x,\theta)-b/a$. Here, the complex coefficients $C_{1}$ and $C_{2}$
represent the amplitudes and phases of the corresponding modes, which are
given by $|C_{i}|$ and
$\Theta_{i}=\tan^{-1}\big{(}\imaginary(C_{i})/\real(C_{i})\big{)}$,
respectively. $\bar{C}$ indicates the complex conjugate of $C$. As
$\bm{u}(t,x,\theta)$ is a real function, the coefficients of the pairwise
eigenfunction are complex conjugates. Notably, a single pair of
eigenfunctions, $W_{\lambda,\pm k}(x,\theta)=X_{\lambda,\pm k}(x)e^{\pm
ik\theta}$, is insufficient to express the helical pattern; at least two pairs
of eigenfunctions are required. In the present example, the two eigenvalues
are indicated by the colored dots in the dispersion relation shown in Fig.
3(d), and the corresponding eigenmodes are shown in Fig. 3(e). As the pattern
propagates along the $\theta$-axis, the coefficients of these modes, $C_{1}$
and $C_{2}$, oscillate periodically with constant amplitudes $|C_{1}|$ and
$|C_{2}|$, respectively, and with the same frequency $k\omega$, where $\omega$
is the angular velocity of the pattern propagation (Figs. 3(f) and (g)). Thus,
Eq. (9) can be expressed as

$\displaystyle\bm{U}_{\ast}(t,x,\theta)$
$\displaystyle=2\big{(}|C_{1}|\bm{A}_{\lambda_{1}}X_{1}\cos(k(\theta-\omega
t)+\theta_{1})$
$\displaystyle\quad+|C_{2}|\bm{A}_{\lambda_{2}}X_{2}\cos(k(\theta-\omega
t)+\theta_{2})\big{)}~{},$ (10)

where $\theta_{1}$ and $\theta_{2}$ denote, respectively, the phases of
$C_{1}$ and $C_{2}$ at $t=0$. This expression indicates that the pattern
propagates with a constant shape and velocity.

The static patterns are also well approximated by two pairwise eigenfunctions
with the same wavenumber $k$, where $C_{1}$ and $C_{2}$ are constants (Figs.
S2–S4 [68]). Overall, these results indicate that the time evolution of the
propagating pattern can be described by the coefficients $C_{i}(t)$ for a few
modes.

###  III.2 Amplitude equations for axisymmetric surfaces

Based on the observations earlier, we approximate a pattern on a surface using
two pairs of eigenfunctions with the same wavenumber $k$ in the
$\theta$-direction as

$\displaystyle\bm{U}(t,x,\theta)$
$\displaystyle=C_{1}\bm{A}_{\lambda_{1}}W_{\lambda_{1},k}+C_{2}\bm{A}_{\lambda_{2}}W_{\lambda_{2},k}+c.c.~{},$
(11)

where $c.c.$ indicates the complex conjugate of the former terms, and the
coefficients $C_{1}(t)$ and $C_{2}(t)$ are functions of time $t$. We explore
amplitude equations that describe the pattern dynamics on an axisymmetric
curved surface.

$\displaystyle\partial_{t}C_{1}$ $\displaystyle=H_{1}(C_{1},C_{2})~{},$ (12)
$\displaystyle\partial_{t}C_{2}$ $\displaystyle=H_{2}(C_{1},C_{2})~{}.$ (13)

We mainly consider the situation in which two pairs of real eigenfunctions
with the same $k$, namely $\bm{A}_{\lambda_{1}}W_{\lambda_{1},k}$ and
$\bm{A}_{\lambda_{2}}W_{\lambda_{1},k}$, are relevant. This situation meets
the aforementioned helical pattern and other cases (see Supplemental Text Sec.
S4 for cases in which the two modes have different wavenumbers $k_{1}$ and
$k_{2}$ [68]). In the vicinity of the Turing bifurcation point, the amplitudes
of these modes, $|C_{1}|$ and $|C_{2}|$, are both expected to be small, and
the right-hand sides of the amplitude equations, $H_{1}(C_{1},C_{2})$ and
$H_{2}(C_{1},C_{2})~{}$, can be expanded with respect to $C_{1}$ and $C_{2}$
[63]. The symmetries inherent in the system impose strong restrictions on the
possible forms of the amplitude equations, as investigated in this section.
The amplitude equations can also be derived from the original RDE using
methods such as the reductive perturbation method [64]; we will carry out the
derivation for the Brusselator model in Sec. V.

####  III.2.1 Amplitude equations for general axisymmetric surfaces

For an axisymmetric surface, the RDE is invariant about the translation by
arbitrary $p_{\theta}$ along the $\theta$-direction, $\mathcal{T}_{\theta}$,
and the reflection about the $\theta$-direction, $\mathcal{M}_{\theta}$:

$\displaystyle\mathcal{T}_{\theta}$
$\displaystyle:\theta\mapsto\theta+p_{\theta},~{}^{\forall}p_{\theta}\in\mathbb{R}~{},$
(14) $\displaystyle\mathcal{M}_{\theta}$
$\displaystyle:\theta\mapsto-\theta~{}.$ (15)

The actions of the transformations $\mathcal{T}_{\theta}$ and
$\mathcal{M}_{\theta}$ on a pattern (see Eq. (11)) can be recast because their
actions on the amplitudes are as follows:

$\displaystyle\mathcal{T}_{\theta}$
$\displaystyle:(C_{1},C_{2})\mapsto(C_{1}e^{-ikp_{\theta}},C_{2}e^{-ikp_{\theta}})~{},$
(16) $\displaystyle\mathcal{M}_{\theta}$
$\displaystyle:(C_{1},C_{2})\mapsto(\bar{C}_{1},\bar{C}_{2})~{}.$ (17)

The amplitude equations, which are the reduction equations of RDE, should also
be invariant to these transformations. Thus, the general form of the amplitude
equations up to the third-order expansion of $C_{i}$ is given by

$\displaystyle\partial_{t}C_{1}$
$\displaystyle=a_{1}C_{1}+a_{2}|C_{1}|^{2}C_{1}+a_{3}|C_{2}|^{2}C_{2}+a_{4}|C_{1}|^{2}C_{2}$
$\displaystyle\quad+a_{5}|C_{2}|^{2}C_{1}+a_{6}C_{1}^{2}\bar{C}_{2}+a_{7}C_{2}^{2}\bar{C}_{1}~{},$
(18) $\displaystyle\partial_{t}C_{2}$
$\displaystyle=b_{1}C_{2}+b_{2}|C_{1}|^{2}C_{1}+b_{3}|C_{2}|^{2}C_{2}+b_{4}|C_{1}|^{2}C_{2}$
$\displaystyle\quad+b_{5}|C_{2}|^{2}C_{1}+b_{6}C_{1}^{2}\bar{C}_{2}+b_{7}C_{2}^{2}\bar{C}_{1}~{},$
(19)

where the coefficients $a_{j}$ and $b_{j}$ are real values owing to the
invariance of the transformation $\mathcal{M}_{\theta}$. The second-order
terms do not appear in the equations. The coefficients $a_{1}$ and $b_{1}$
coincide with growth rates $\mu(\lambda_{1})$ and $\mu(\lambda_{2})$,
respectively. The other coefficients are not determined solely by the
discussion of symmetries given here. Below, we assume that the coefficients
are chosen such that the solution does not diverge; otherwise, higher-order
terms are required to prevent divergence.

Using the complex variables $C_{1}$ and $C_{2}$ as $C_{1}=\eta
e^{i\phi},~{}C_{2}=\xi e^{i\psi}$, the pattern profile in Eq. (11) can be
written as

$\displaystyle\bm{U}(t,x,\theta)$
$\displaystyle=2\eta\bm{A}_{\lambda_{1}}X_{\lambda_{1},k}\cos(k\theta+\phi)$
$\displaystyle\quad+2\xi\bm{A}_{\lambda_{2}}X_{\lambda_{2},k}\cos(k\theta+\psi)~{},$
(20)

where the amplitudes, $\eta$ and $\xi$, and phase variables, $\phi$ and
$\psi$, are real-valued functions of time $t$. The amplitude equations can be
rewritten as follows:

$\displaystyle\partial_{t}\eta$
$\displaystyle=\big{(}a_{1}+a_{2}\eta^{2}+a_{5}\xi^{2}\big{)}\eta+a_{7}\eta\xi^{2}\cos(2\alpha)$
$\displaystyle\quad+\big{(}a_{3}\xi^{2}+(a_{4}+a_{6})\eta^{2}\big{)}\xi\cos(\alpha)~{},$
(21) $\displaystyle\partial_{t}\xi$
$\displaystyle=\big{(}b_{1}+b_{3}\xi^{2}+b_{4}\eta^{2}\big{)}\xi+b_{6}\eta^{2}\xi\cos(2\alpha)$
$\displaystyle\quad+\big{(}b_{2}\eta^{2}+(b_{5}+b_{7})\xi^{2}\big{)}\eta\cos(\alpha)~{},$
(22) $\displaystyle\eta\xi\partial_{t}\alpha$
$\displaystyle=-\big{(}a_{3}\xi^{4}+(a_{4}-a_{6}+b_{5}-b_{7})\eta^{2}\xi^{2}+b_{2}\eta^{4}$
$\displaystyle\qquad\quad+2(a_{7}\xi^{2}+b_{6}\eta^{2})\eta\xi\cos(\alpha)\big{)}\sin(\alpha)~{},$
(23)

where we introduce the phase difference $\alpha\equiv\phi-\psi$. These
equations are closed in three variables, $\eta$, $\xi$, and $\alpha$, because
the system is invariant for translation along the $\theta$-axis. Note also
that they are invariant under the transformation $\alpha\rightarrow-\alpha$
because the system has reflection-symmetry along the $\theta$-axis. The phase
variables $\phi$ and $\psi$ evolve according to the following equations:

$\displaystyle\eta\partial_{t}\phi$
$\displaystyle=-\xi((a_{4}-a_{6})\eta^{2}+a_{3}\xi^{2})\sin(\alpha)$
$\displaystyle\quad-a_{7}\eta\xi^{2}\sin(2\alpha),$ (24)
$\displaystyle\xi\partial_{t}\psi$
$\displaystyle=+\eta(b_{2}\eta^{2}+(b_{5}-b_{7})\xi^{2})\sin(\alpha)$
$\displaystyle\quad+b_{6}\eta^{2}\xi\sin(2\alpha).$ (25)

Eqs. (21)–(23), supplemented by Eqs. (24) and (25), constitute the general
form of the amplitude equations for pattern dynamics on an axisymmetric
surface. Note that the derived amplitude equations Eqs. (18) and (19) are also
valid for non-Turing systems as long as the system satisfies similar
conditions such as invariance to the transformations Eqs. (14) and (15).

####  III.2.2 Amplitude equations for reflection-symmetric surfaces

An axisymmetric surface can have additional symmetries. Here, we consider the
case in which the surface is reflection-symmetric about the $x$-axis,
satisfying $r(x)=r(-x)$, where the origin of the $x$-axis can be arbitrary
because of the periodic boundary condition. For such a surface, RDE is
invariant under the reflection about the $x$-direction, $\mathcal{M}_{x}$:

$\displaystyle\mathcal{M}_{x}$ $\displaystyle:x\mapsto-x~{}.$ (26)

The eigenfunctions of the Laplace–Beltrami operator $X_{\lambda,k}(x)$ are
either even or odd with respect to $x$. The action of the reflection
$\mathcal{M}_{x}$ on $C_{i}$ is dependent on the parity of
$X_{\lambda_{i},k}(x)$ and acts as

$\displaystyle\mathcal{M}_{x}$
$\displaystyle:C_{i}\mapsto\left\\{\begin{array}[]{cl}+C_{i}&~{}\mbox{if
$X_{\lambda_{i},k}$ is an even function}~{},\\\ -C_{i}&~{}\mbox{if
$X_{\lambda_{i},k}$ is an odd function}~{}.\end{array}\right.$ (29)

When $X_{\lambda_{1},k}(x)$ and $X_{\lambda_{2},k}(x)$ are both even or odd
functions, the amplitude equations are the same as those in the general case
of Eqs. (18) and (19). However, when $X_{\lambda_{1},k}(x)$ and
$X_{\lambda_{2},k}(x)$ are a combination of even and odd functions, the
amplitude equations are reduced to

$\displaystyle\partial_{t}C_{1}$
$\displaystyle=a_{1}C_{1}+a_{2}|C_{1}|^{2}C_{1}+a_{5}|C_{2}|^{2}C_{1}+a_{7}C_{2}^{2}\bar{C}_{1}~{},$
(30) $\displaystyle\partial_{t}C_{2}$
$\displaystyle=b_{1}C_{2}+b_{3}|C_{2}|^{2}C_{2}+b_{4}|C_{1}|^{2}C_{2}+b_{6}C_{1}^{2}\bar{C}_{2}~{},$
(31)

corresponding to the case $a_{3}=a_{4}=a_{6}=0$ and $b_{2}=b_{5}=b_{7}=0$ in
Eqs. (18) and (19). The amplitude and phase variables obey

$\displaystyle\partial_{t}\eta$
$\displaystyle=\big{(}a_{1}+a_{2}\eta^{2}+a_{5}\xi^{2}\big{)}\eta+a_{7}\eta\xi^{2}\cos(2\alpha)~{},$
(32) $\displaystyle\partial_{t}\xi$
$\displaystyle=\big{(}b_{1}+b_{3}\xi^{2}+b_{4}\eta^{2}\big{)}\xi+b_{6}\eta^{2}\xi\cos(2\alpha)~{},$
(33) $\displaystyle\partial_{t}\alpha$
$\displaystyle=-(a_{7}\xi^{2}+b_{6}\eta^{2})\sin(2\alpha)~{},$ (34)

and

$\displaystyle\partial_{t}\phi$ $\displaystyle=-a_{7}\xi^{2}\sin(2\alpha)~{},$
(35) $\displaystyle\partial_{t}\psi$
$\displaystyle=+b_{6}\eta^{2}\sin(2\alpha)~{}.$ (36)

####  III.2.3 Amplitude equations for $2\pi$-periodic surface along the
$x$-axis

The axisymmetric surface can be periodic along the $x$-axis, satisfying
$r(x)=r(x+p_{x})$, where $p_{x}$ is a constant smaller than the axial length
of the surface. Here, we consider a simple case in which the surface is
$2\pi$-periodic (as in our previous study [57]), satisfying $r(x)=r(x+2\pi)$,
and the domain of $x$ is defined as $-2\pi\leq x<2\pi$. For such a surface,
the RDE is invariant to the $2\pi$-translation along the $x$-axis and
$\mathcal{T}_{x,2\pi}$ is defined as

$\displaystyle\mathcal{T}_{x,2\pi}$ $\displaystyle:x\mapsto x+2\pi~{},$ (37)

and reflection about the $\theta$-direction, $\mathcal{M}_{\theta}$. According
to Bloch’s theorem [71], the eigenfunctions of the Laplace–Beltrami operator
$X_{\lambda,k}(x)$ either change sign or remain unchanged by the action of
$\mathcal{T}_{x,2\pi}$.

$\displaystyle X_{\lambda,k}(x+2\pi)=\pm X_{\lambda,k}(x).$ (38)

The action of $2\pi$-translation $\mathcal{T}_{x,2\pi}$ on $C_{i}$ act as

$\displaystyle\mathcal{T}_{x,2\pi}$
$\displaystyle:C_{i}\mapsto\left\\{\begin{array}[]{cl}+C_{i}&~{}\mbox{if
$X_{\lambda_{i},k}(x+2\pi)=+X_{\lambda_{i},k}(x)$},\\\ -C_{i}&~{}\mbox{if
$X_{\lambda_{i},k}(x+2\pi)=-X_{\lambda_{i},k}(x)$}\end{array}\right.$ (41)

When $X_{\lambda_{1},k}(x)$ and $X_{\lambda_{2},k}(x)$ both undergo the same
transformation by $\mathcal{T}_{x,2\pi}$, the amplitude equations are the same
as those in the general case of Eqs. (18) and (19). Conversely, when
$X_{\lambda_{1},k}(x)$ and $X_{\lambda_{2},k}(x)$ undergo transformations of
different signs by $\mathcal{T}_{x,2\pi}$, the amplitude equations are the
same as those on the reflection-symmetric surface in Eqs. (30) and (31).

####  III.2.4 Patterns consisting of single paired modes

The patterns can be dominated by only a single paired mode, as exemplified in
Fig. 4, where the amplitudes of the other modes decrease.

$\displaystyle\bm{U}=C\bm{A}_{\lambda}X_{\lambda,k}e^{ik\theta}+c.c.~{}.$ (42)

The pattern is always reflection-symmetric around the $\theta$-axis as
$U(x,\theta)=U(x,-\theta)$ where the origin of the $\theta$-axis is properly
chosen. Based on the invariance by the transformations $\mathcal{T}_{\theta}$
and $\mathcal{M}_{\theta}$, the amplitude equation is given by

$\displaystyle\partial_{t}C$ $\displaystyle=a_{1}C+a_{2}|C|^{2}C~{},$ (43)

where $a_{1}$ coincides with the growth rate $\mu(\lambda)$ and $a_{2}$ is
supposed to be a negative real value. With $C=\eta e^{i\phi}$, we obtain:

$\displaystyle\partial_{t}\eta$ $\displaystyle=a_{1}\eta+a_{2}\eta^{3}~{},$
(44) $\displaystyle\partial_{t}\phi$ $\displaystyle=0~{}.$ (45)

These equations indicate that the pattern does not exhibit propagation.

####  III.2.5 Patterns on a flat surface

Summarizing the analysis for the flat surface case is helpful. On the
$xy$-plane, the pattern can be expanded by using the Fourier mode
$e^{i(k_{x}x+k_{y}y)}$. The stripe pattern resulting from Turing instability
(Fig. 2) comprises at most two paired modes and is approximated by

$\displaystyle\bm{U}=(C_{1}e^{i(k_{x}x+k_{y}y)}+C_{2}e^{i(k_{x}x-k_{y}y)})\bm{A}+c.c.~{}.$
(46)

The amplitude equations are derived by considering that the plane is invariant
to translations and reflections:

$\displaystyle\partial_{t}C_{1}$
$\displaystyle=a_{1}C_{1}+a_{2}|C_{1}|^{2}C_{1}+a_{5}|C_{2}|^{2}C_{1}~{},$
(47) $\displaystyle\partial_{t}C_{2}$
$\displaystyle=b_{1}C_{2}+b_{3}|C_{2}|^{2}C_{2}+b_{4}|C_{1}|^{2}C_{2}~{}.$
(48)

Then, amplitudes and phases obey

$\displaystyle\partial_{t}\eta$
$\displaystyle=a_{1}\eta+a_{2}\eta^{3}+a_{5}\eta\xi^{2}~{},$ (49)
$\displaystyle\partial_{t}\xi$
$\displaystyle=b_{1}\xi+b_{2}\xi^{3}+b_{4}\eta^{2}\xi~{},$ (50)
$\displaystyle\partial_{t}\phi$ $\displaystyle=\partial_{t}\psi=0~{},$ (51)

and pattern propagation is not possible. In the case of a plane, coupling with
higher-order harmonic modes such as $e^{ik_{x}x}$ and $e^{i2k_{x}x}$ can
result in propagating waves, even in a one-dimensional (1D) system [63]. We
ignore such situations because they generally occur far from the Turing
bifurcation point and depend on the specific nature of the system, such as
domain size (see Supplemental Text Sec. S4 [68]).

##  IV Steady-state solutions of amplitude equations

Here, we investigate the steady-state solutions of the amplitude equations
derived in the preceding section:
$(\eta,\xi,\alpha)=(\eta_{0},\xi_{0},\alpha_{0})$. For the given solution, the
symmetries of the corresponding patterns and phase velocities
$\partial_{t}\phi,\partial_{t}\psi$ are determined. As a pattern of single
paired modes cannot propagate (see Eq. (45)), we consider the case in which
two pairs of modes are relevant and $\eta_{0}$ and $\xi_{0}$ are assumed to be
non-zero.

###  IV.1 Steady-state solutions of Eqs. (32)–(34)

The amplitude equations (32)–(34) are obtained for the reflection-symmetric
surface with eigenfunctions that have different parities as
$X_{1}(-x)=X_{1}(x)$ and $X_{2}(-x)=-X_{2}(x)$ (Sec. III.2.2), or for a
$2\pi$-periodic surface with eigenfunctions that transform differently as
$X_{1}(x+2\pi)=X_{1}(x)$ and $X_{2}(x+2\pi)=-X_{2}(x)$ (Sec. III.2.3). In the
steady-state, either of the following relationships should hold from Eq. (34):

$\displaystyle\alpha_{0}=\frac{n\pi}{2}~{},~{}n\in\mathbb{Z}~{},$ (52)
$\displaystyle a_{7}\xi_{0}^{2}+b_{6}\eta_{0}^{2}=0~{}.$ (53)

For the solution $\alpha_{0}=n\pi/2$ (Eq. (52)), the phase equations in Eqs.
(35) and (36) satisfy $\partial_{t}\phi=\partial_{t}\psi=0$, indicating that
this solution represents a static pattern. The profile of the pattern is
expressed as

$\displaystyle{\bm{U}}(x,\hat{\theta})$
$\displaystyle=2\eta_{0}\bm{A}_{\lambda_{1}}X_{\lambda_{1},k}(x)\cos\bigl(\hat{\theta}\bigr{missing})+2\xi_{0}\bm{A}_{\lambda_{2}}X_{\lambda_{2},k}(x)\cos\bigl(\hat{\theta}-\alpha_{0}\bigr{missing})~{}$
(54)
$\displaystyle=\begin{dcases}2\eta_{0}\bm{A}_{\lambda_{1}}X_{\lambda_{1},k}(x)\cos\bigl(\hat{\theta}\bigr{missing})+2(-1)^{n-1}\xi_{0}\bm{A}_{\lambda_{2}}X_{\lambda_{2},k}(x)\cos\bigl(\hat{\theta}\bigr{missing})&\text{if
$\alpha_{0}=n\pi$}\\\
2\eta_{0}\bm{A}_{\lambda_{1}}X_{\lambda_{1},k}(x)\cos\bigl(\hat{\theta}\bigr{missing})+2(-1)^{n-1}\xi_{0}\bm{A}_{\lambda_{2}}X_{\lambda_{2},k}(x)\sin\bigl(\hat{\theta}\bigr{missing})&\text{if
$\alpha_{0}=(n+1/2)\pi$}\end{dcases}$ (55)

where $\hat{\theta}=k\theta+\phi_{0}$ and $\phi_{0}$ is a constant determined
by the initial condition. When $\alpha_{0}=n\pi$, the pattern is reflection-
symmetric along the $\theta$-axis, satisfying
${\bm{U}}(x,\hat{\theta})={\bm{U}}(x,-\hat{\theta})$. Fig. 5(a) shows an
example of a $2\pi$-periodic surface, where we set
$r(x)=2+0.3\cos(x)+0.2\sin(x)$. In the other case, $\alpha_{0}=(n+1/2)\pi$,
the pattern symmetry depends on the surface symmetry. On reflection-symmetric
surfaces, the pattern is inversion-symmetric, satisfying
${\bm{U}}(x,\hat{\theta})={\bm{U}}(-x,-\hat{\theta})$. On $2\pi$-periodic
surfaces, the pattern is invariant against the simultaneous transformation of
reflection about the $\theta$-axis and $2\pi$-translation along the $x$-axis,
satisfying ${\bm{U}}(x,\hat{\theta})={\bm{U}}(x+2\pi,-\hat{\theta})$. Fig.
5(b) shows an example of the case on a reflection-symmetric surface
($r(x)=2.2+0.2\cos(x)+0.4\cos(x/2)$).

For the steady-state solution given by Eq. (53), the phase velocities
$\partial_{t}\phi$ and $\partial_{t}\psi$ are generally non-zero and constant
($\partial_{t}\phi=\partial_{t}\psi\neq 0$), indicating a propagating pattern.
The pattern shows neither reflection- nor inversion- symmetry because of the
superposition of the two paired modes in $\alpha_{0}\neq n\pi/2$. This non-
trivial solution was omitted in the previous approach [57]. Fig. 5(c) shows an
example of a propagating pattern on a $2\pi$-periodic surface
($r(x)=1.9+0.21\cos(x)+0.06\sin(2x)$). This non-trivial solution is absent on
a flat-plane where $b_{6}=a_{7}=0$ (see Eqs. (47) and (48)).

###  IV.2 Steady-state solutions of Eqs. (21)–(23)

Next, we consider the steady-state solutions of the general amplitude
equations given by Eqs. (21)–(23). The equations are obtained for a general
axisymmetric surface (Sec. III.2.1), for cases in which $X_{1}(x)$ and
$X_{2}(x)$ are both even or odd functions on a reflection-symmetric surface
(Sec. III.2.2), or when $X_{1}(x)$ and $X_{2}(x)$ are transformed with the
same sign by $2\pi$-translation along the $x$-axis on a $2\pi$-periodic
surface (Sec. III.2.3). In the steady-state, Eq. (23) yields to either of the
following solutions:

$\displaystyle\alpha_{0}=n\pi,~{}n\in\mathbb{Z},$ (56)
$\displaystyle\cos\alpha_{0}=-\frac{a_{3}\xi_{0}^{4}+(a_{4}-a_{6}+b_{5}-b_{7})\eta_{0}^{2}\xi_{0}^{2}+b_{2}\eta_{0}^{4}}{2(a_{7}\xi_{0}^{2}+b_{6}\eta_{0}^{2})\eta_{0}\xi_{0}}~{},$
(57) $\displaystyle
a_{3}\xi_{0}^{4}+(a_{4}-a_{6}+b_{5}-b_{7})\eta_{0}^{2}\xi_{0}^{2}+b_{2}\eta_{0}^{4}=0$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mbox{with}~{}~{}a_{7}\xi_{0}^{2}+b_{6}\eta_{0}^{2}=0~{}.$
(58)

From Eqs. (24) and (25), the first solution $\alpha_{0}=n\pi$ represents a
static pattern whereas the others represent propagating patterns. Note that at
$a_{3}=a_{4}=a_{6}=b_{2}=b_{5}=b_{7}=0$, the solutions Eqs. (57) and (58)
coincide with the solutions $\alpha_{0}=(n+1/2)\pi$ included in Eq. (52) and
Eq. (53), respectively.

For the solution $\alpha_{0}=n\pi$, the static pattern takes the form of the
first line of Eq. (55) and is reflection-symmetric about the $\theta$-axis.
The pattern has an additional symmetry depending on the surface symmetry and
combination of eigenfunctions. When both $X_{1}(x)$ and $X_{2}(x)$ are even
functions on the reflection-symmetric surface, the pattern has symmetry
${\bm{U}}(x,\hat{\theta})={\bm{U}}(x,-\hat{\theta})={\bm{U}}(-x,\hat{\theta})={\bm{U}}(-x,-\hat{\theta})$,
while when both are odd, it has the symmetry
${\bm{U}}(x,\hat{\theta})={\bm{U}}(x,-\hat{\theta})={\bm{U}}(-x,\hat{\theta}+\pi)={\bm{U}}(-x,-\hat{\theta}+\pi)$.
An example of the latter static pattern is shown in Fig. 6(a)
($r(x)=1.95+0.21\cos(x)-0.06\cos(2x)+0.015\cos(x/2)$). On a $2\pi$-periodic
surface, when $X_{1}(x+2\pi)=X_{1}(x)$ and $X_{2}(x+2\pi)=X_{2}(x+2\pi)$, the
pattern has the symmetry
${\bm{U}}(x,\hat{\theta})={\bm{U}}(x,-\hat{\theta})={\bm{U}}(x+2\pi,\hat{\theta})={\bm{U}}(x+2\pi,-\hat{\theta})$,
while when $X_{1}(x+2\pi)=-X_{1}(x)$ and $X_{2}(x+2\pi)=-X_{2}(x+2\pi)$, it
has
${\bm{U}}(x,\hat{\theta})={\bm{U}}(x,-\hat{\theta})={\bm{U}}(x+2\pi,\hat{\theta}+\pi)={\bm{U}}(x+2\pi,-\hat{\theta}+\pi)$.
An example of the latter static pattern is shown in Fig. 6(b)
($r(x)=1.8+0.15\cos(x)+0.09\cos(2x)$).

For other solutions, Eqs. (57) and (58), $\alpha_{0}$ takes a non-trivial
value, and the corresponding propagating patterns have no symmetry. Fig. 6(c)
presents an example for the solution of Eq. (57) on a general axisymmetric
surface ($r(x)=2.2+0.2\cos(x)+0.4\cos(x/2-\pi/4)$). Regarding the solution Eq.
(58), the parameters are required to satisfy additional conditions, such as
$a_{3}b_{6}^{2}+b_{2}a_{7}^{2}-a_{7}b_{6}(a_{4}-a_{6}+b_{5}-b_{7})=0$, which
are unlikely to be achieved.

###  IV.3 Classification of parallel and helical patterns by steady-state
solution of amplitude equations

The dynamics of the parallel and helical patterns (Fig. 1) can be understood
using the steady-state solutions discussed above. On a $2\pi$-periodic surface
along the $x$-axis, parallel patterns are static, whereas helical patterns are
static only on a reflection-symmetric surface (Fig. 7 (a)–(c)). These patterns
correspond to the steady-state solution $\alpha_{0}=\pi/2$ in Eq. (52), each
consisting of eigenfunctions with characteristic symmetry; for instance, the
helical pattern on a reflection-symmetric surface is composed of
eigenfunctions with $X_{1}(x)=X_{1}(-x)=X_{1}(x+2\pi)$ and
$X_{2}(x)=-X_{2}(-x)=X_{2}(x+2\pi)$. The helical pattern begins to move
because of the loss of reflection-symmetry [57], which corresponds to the
steady-state solution in Eq. (57) of the amplitude equations in Eqs.
(21)–(23), where $X_{1}(x)=X_{1}(x+2\pi)$, $X_{2}(x)=-X_{2}(x+2\pi)$ and the
phase difference $\alpha_{0}\neq n\pi$ (Fig. 7(d)).

The parallel static pattern on a $2\pi$-periodic surface along the $x$-axis
satisfies $X_{1}(x)=X_{1}(x+2\pi)$ and $X_{2}(x)=-X_{2}(x+2\pi)$. Because of
the loss of the $2\pi$-periodicity of the surface, these relationships no
longer hold, and the solution of the amplitude equations becomes Eq. (52),
indicating the onset of a propagating wave in a parallel pattern. A numerical
example is shown in Fig. 8.

###  IV.4 Summary of amplitude equations and steady-state solutions

We derived the amplitude equations (18) and (19) for a general axisymmetric
surface on which two pairs of modes couple with each other, and the static and
propagating patterns correspond to the steady solutions of the amplitude
equations. Some coupling terms are eliminated because of the symmetry of the
surface (see Eqs. (30) and (31)), which suppresses pattern propagation. The
propagation of helical and parallel patterns is caused by the loss of
reflection-symmetry and $2\pi$-periodicity, respectively, of the surface along
the $x$-axis. For a flat surface, the amplitude equations are reduced to Eqs.
(47) and (48), and the pattern is always static. The relationships obtained
for the pattern dynamics and the symmetry of the surface and pattern are
consistent with those reported in a previous study [57]. In addition, from the
amplitude equations, we obtained a non-trivial solution (Eq. (53)) for which
the pattern generally has no specific symmetry and propagates. This solution
allows the Turing pattern to propagate even on highly symmetric surfaces. A
propagating pattern on a sphere was observed in numerical simulation [57].

##  V Amplitude equation for the Brusselator model

###  V.1 Reductive perturbation method

The amplitude equations for RDE on axisymmetric curved surfaces can be derived
by applying the reductive perturbation method to reaction-diffusion systems in
the vicinity of the Turing bifurcation. Consider the general RDE expressed in
Eq. (1), for which we explicitly include a bifurcation parameter $\nu$ in the
reaction term as $\bm{R}({\bm{u}},\nu)$. As before, we assume that the
equation has a unique uniform steady solution, $\bm{u}_{0}$, satisfying
$\bm{R}(\bm{u}_{0},\nu)=0$. The uniform solution becomes unstable for $\nu>0$,
and a Turing pattern occurs.

In the vicinity of the Turing bifurcation at a positive, sufficiently small
value of $\nu$, the amplitude of the pattern,
$|{\bm{U}}|\equiv|{\bm{u}}-{\bm{u}}_{0}|$, is expected to be small, and
expanding the RDE by $\bm{U}$ and $\nu$ is feasible. The amplitude equations
can then be derived from the RDE; see the details of the derivation in
Supplemental Text Sec. S2 [68]. In summary, the pattern dynamics of the RDE
are approximated by a lower-order expansion with respect to $\nu$ and
amplitude $C_{i}$, and the pattern dynamics are dominated by the time
evolution of $C_{i}$:

$\displaystyle\bm{U}(t)$
$\displaystyle=\bm{U}_{0}(\\{C_{i}(t)\\})+\bm{\rho}(\\{C_{i}(t)\\})~{},$ (59)
$\displaystyle\partial_{t}C_{i}(t)$ $\displaystyle=H_{i}(\\{C_{i}(t)\\})~{},$
(60)

where $\bm{U}_{0}$ and $\bm{\rho}$ are the first- and higher-order terms in
the amplitudes, respectively, and the second equation is referred to as the
amplitude equation. These are derived in Eq. (S41) and in Eq. (S42),
respectively. The obtained equation (S42) not only validates Eqs. (18) and
(19), which are derived as a general form of the amplitude equations, but also
determines their coefficients. The simpler forms of amplitude equations (30)
and (31) can be understood as follows. As shown in Eqs. (S44)–(S49), the
coefficients include integrals of the products of the Laplace–Beltrami
eigenfunctions $X_{\lambda_{i},k}$. For example, if the surface is reflection-
symmetric about the $x$-axis and $X_{1}(x)$ and $X_{2}(x)$ are a pair of even
and odd functions, some coefficients are equal to zero owing to the
combination of the parity of the eigenfunctions in the integrands, and Eq.
(S42) coincides with Eqs. (30) and (31).

###  V.2 Derivation of Amplitude equation for Brusselator model

We derive the amplitude equations for the Brusselator model in Eqs. (18) and
(19) following the general procedure described in Supplemental Text Sec. S2
[68]. We choose the bifurcation parameter $\nu=b-b_{0}$, where the reference
parameters of the reductive perturbation, $(a,b_{0},D_{u},D_{v})$, are set
such that the growth rates at the two Laplace–Beltrami eigenvalues,
$\mu(\lambda_{1})$ and $\mu(\lambda_{2})$, both become zero. For positive
$\nu$, the pattern can be approximated as

$\displaystyle{\bm{U}}=C_{1}(t){\bm{A}}_{\lambda_{1}}W_{\lambda_{1},k}+C_{2}(t){\bm{A}}_{\lambda_{2}}W_{\lambda_{2},k}+c.c.~{},$
(61)

where ${\bm{A}}_{\lambda_{1}}$ and ${\bm{A}}_{\lambda_{2}}$ are eigenvectors
of the linear operator:

$\displaystyle\mathcal{L}_{0}=\begin{pmatrix}-\lambda_{i}D_{u}+b_{0}-1&a^{2}\\\
-b_{0}&-\lambda_{i}D_{v}-a^{2}\end{pmatrix}~{}.$ (62)

Here, $\lambda_{i}~{}(i=1,2)$ are the eigenvalues of the Laplace–Beltrami
operator for $W_{\lambda_{i},k}$ and are determined from the surface shape.
The eigenvalues $\lambda_{i}$ and parameters $(a,b_{0},D_{u},D_{v})$ satisfy
the relation explained in Supplemental Text Sec. S3 [68]. The eigenvectors
${\bm{A}}_{\lambda_{i}}$ and their conjugate vectors ${\bm{B}}_{\lambda_{i}}$
are expressed by Eqs. (S52) and (S53), respectively. The coefficients of the
amplitude equations are then calculated using the products of
$\bm{A}_{\lambda_{i}}$ and $\bm{B}_{\lambda_{i}}$, and the integral with
respect to $X_{\lambda_{i},k}$. For more details, see Supplemental Text Sec.
S3 [68].

###  V.3 Comparison between Brusselator model and amplitude equations

We numerically solve the amplitude equations derived earlier, and we compare
the solution with that of the Brusselator model. In particular, we explore the
steady-state solution of the amplitude equations, which corresponds to either
a static or propagating pattern, and confirm the claims given in Sec. IV.
Dynamics other than the static and propagating patterns are discussed in the
next section.

####  V.3.1 Comparison procedure

The amplitude equations were derived using the following procedure. First, we
fixed the surface and selected two pairs of Laplace–Beltrami eigenvalues,
$\lambda_{1}$ and $\lambda_{2}$, that were adjacent in the spectral space, and
then determined the corresponding eigenfunctions, $W_{1,\pm k}$ and $W_{2,\pm
k}$. These were chosen to have the same wavenumber $k$. Next, we set the
diffusion coefficients $D_{u}$ and $D_{v}$ and then determined the model
parameters $a$ and $b_{0}$ to satisfy the determinant of $\mathcal{L}_{0}$ to
zero (see Eqs. (S58) and (S59)). Subsequently, we set the bifurcation
parameter $\nu(=b-b_{0})$, which was chosen to be sufficiently small and of
the order of $10^{-2}$. The parameters were chosen so that the system would
exhibit Turing instability, where the growth rates $\mu(\lambda_{1})$ and
$\mu(\lambda_{2})$ were real and positive values. Finally, we determined the
coefficients of the amplitude equations, $a_{j}$ and $b_{j}$
($j=1,2,\ldots,7$), from Eqs. (S43)–(S49), and performed numerical
simulations. We also numerically solved the Brusselator model (the original
system) using the same surface and parameters.

We compared the solutions obtained from the amplitude equations and the
Brusselater model in the phase space of $\eta$, $\xi$ and phase differences
$\alpha$. For the Brusselator model, $(\eta,\xi,\alpha)$ was obtained by
taking the inner product between the pattern $\bm{U}$ and eigenvectors
$\bm{B}_{\lambda_{1}}W_{\lambda_{1},\pm k}$ and
$\bm{B}_{\lambda_{2}}W_{\lambda_{2},\pm k}$, where
$\bm{U}(t,x,\theta)=\bm{U}_{0}(t,x,\theta)+\mathcal{O}(|\nu|)$ and

$\displaystyle\bm{U}_{0}(t,x,\theta)$
$\displaystyle=\eta(t)e^{i\phi(t)}\bm{A}_{\lambda_{1}}W_{\lambda_{1},k}(x,\theta)$
$\displaystyle\quad+\xi(t)e^{i\psi(t)}\bm{A}_{\lambda_{2}}W_{\lambda_{2},k}(x,\theta)+c.c.~{},$
(63)

where $\bm{A}_{\lambda_{1}}$ and $\bm{A}_{\lambda_{2}}$ are the same as the
eigenvectors in the amplitude equations.

####  V.3.2 Comparison of static and propagating patterns

The results of the numerical simulations for the amplitude equations and
Brusselator model are summarized in Fig. 9, in which the blue and orange
points represent the steady solutions $(\eta_{0},\xi_{0},\alpha_{0})$ of the
amplitude equations and Bursselator model, respectively. Patterns on
axisymmetric surfaces and their projection onto the $x$-$\theta$ plane are
also shown, where the left panels show the amplitude equations and the right
panels show the Brusselator model. The coefficients of the amplitude equations
are summarized in Supplemental Table S1 [68].

Figs. 9(a) and (b) show static patterns of the steady solutions. For Fig.
9(a), helical static patterns on the reflection-symmetric and $2\pi$-periodic
surface are obtained through the direct simulation of the Brusselator model,
while the constructed pattern from the amplitude equations exhibits a similar
static pattern with phase difference $\alpha_{0}=\pi/2$. For Fig. 9(b), which
employs the same parameter set as that in Fig. 6(a), the orbit of the
amplitude equations reach a solution similar to that of the corresponding
Brusselator model, with phase difference $\alpha_{0}=0$.

Figs. 9(c) and (d) show propagating patterns with phase differences
$\alpha_{0}\neq n\pi/2$, respectively. For Fig. 9(c), the helical propagating
pattern is constructed from the amplitude equations for a surface without
reflection-symmetry. $\alpha_{0}$ slightly deviates from $\pi/2$ and the
pattern is propagating (see also Fig. 10(c)). The patterns and dynamics
obtained from the amplitude equations and Brusselator model are in close
agreement. For Fig. 9(d), the parameters in Fig. 5(c) are used, where the
phase difference is given by the non-trivial solution Eq. (53). We found that
pattern dynamics constructed by the amplitude equations are similar to those
of the original Brusselator model, indicating that the amplitude equations can
capture the non-trivial propagating solution.

####  V.3.3 Approximation accuracy against bifurcation parameter

The approximation accuracy against the bifurcation parameters of the helical
propagating pattern was evaluated. Under the same conditions as those shown in
Fig. 9(c), the bifurcation parameter $\nu$ was changed from $0.01$ to $0.1$ in
steps of $0.01$, resulting in the observation of similar helical patterns. The
patterns obtained from the amplitude equations closely approximated those of
the Brusselator model, as shown in Figs. 10(a), (b), and (c), in which the
amplitudes $\eta$ and $\xi$ and the phase difference $\alpha$ of amplitude
equations and the original Brusselator model are in close agreement. Figs.
10(a) and (b) confirm that the leading order of the amplitudes is
$\mathcal{O}(|\nu|^{3/2})$, which is assumed to derive the amplitude equations
(Supplemental Text Sec. S2 [68]). The velocities obtained from the amplitude
equations and the original Brusselator model are in close agreement (with
difference $\mathcal{O}(\nu^{3})$), and converge as the bifurcation parameter
$\nu$ approaches zero (Fig. 10(d)). Consequently, the numerical simulations
confirm that the amplitude equations provide a reasonable approximation of the
Brusselator model.

##  VI Pattern dynamics other than propagation

In Sec. III and Sec. V, we derived the amplitude equations for an axisymmetric
surface and investigated the steady-state solutions corresponding to the
static and propagating patterns. The generic form of the derived equations
(Eqs. (18) and (19)) allows us to further explore the possible dynamics, other
than propagation, of the RDE on general axisymmetric surfaces; the pattern
dynamics are explored by numerically searching for non-stationary solutions of
the amplitude equations. In the joint paper [65], we discuss this issue in
detail; this section provides a summary of the results. As explained below,
our results indicate that the Turing pattern can exhibit oscillatory and
chaotic behaviors on curved surfaces.

###  VI.1 Limit cycle solution and pattern dynamics

We explored the possible pattern dynamics on curved surfaces by conducting
numerical searches for amplitude equations (18) and (19). By changing the
surface and model parameters of the Brusselator model, the corresponding
amplitude equations were derived and numerically solved. For a given parameter
set, the derivation of the amplitude equations is the same as that described
in Sec. V.3, where we used the surface
$r(x)=d+\sum_{n=1}^{N}\left(s_{n}\sin(nx)+c_{n}\cos(nx)\right)$ with
$x\in[0,2\pi)$ and $N=5$. For an efficient search, the range of the $x$-axis
was narrowed.

The surface shape and dispersion relations with parameters
$(a,b_{0},D_{u},D_{v})=(1.5,2.97,0.5,2.2)$, $\nu=0.02$,
$d=1,s_{2}=0.15,s_{3}=0.4$ and other surface parameters set to zero are shown
in Fig. 11(a) and (b), respectively. We found that the amplitude equations
have a limit cycle solution, as represented by the blue line in Fig. 11(c)
(see Supplemental Table S1 for the coefficients of the amplitude equations
[68]). The corresponding pattern exhibits propagation accompanied by
oscillations, as shown in Fig. 11(d) left column. The bifurcation analysis
suggests that the oscillatory pattern dynamics appear via Hopf bifurcation
from the propagating pattern [65]. We also validated that the original
Brusselator model yields similar dynamics (Fig.11(c) orange line and (d) right
column). On a flat surface, the pattern remains static as an ordinary Turing
pattern. This example demonstrates that the Turing pattern can exhibit more
complex dynamics than propagation depending on the surface geometry.

###  VI.2 Chaotic dynamics in amplitude equations

We also searched for complex dynamics by changing the coefficients $a_{j}$ and
$b_{j}$ ($j=2,3,\ldots,7$) as arbitrarily selectable parameters, with $a_{1}$
and $b_{1}$ chosen as real positive values to satisfy the condition of Turing
instability. Note that, although this method does not identify the
corresponding Brusselator model because of the difficulty in determining
surface and model parameters from the coefficients, it allows for a general
and efficient investigation of what pattern dynamics are possible. In addition
to the limit cycle solution, we observed chaotic dynamics, as shown in Fig.
12. The appearance of the chaotic dynamics can be understood in terms of
period-doubling bifurcation against parameter $b_{3}$ (Fig. 12(b)). The other
coefficients of the amplitude equations are listed in Supplemental Table S1
[68]. The existence of such a chaotic solution in the amplitude equations
suggests that rich and complex dynamics are possible depending on the geometry
of the surface shape. However, the corresponding patterns have not yet been
obtained for the Brusselator model on axisymmetric surfaces.

As the search using the amplitude equations is limited to the region near the
Turing bifurcation point on axisymmetric surfaces, we examined the Brusselator
model under broader conditions. We found an example of chaotic pattern
dynamics on a deformed spherical surface, as shown in Fig. 13(a) and [65]. The
surface is not axisymmetric and given by a radial function
$R(\theta,\phi)=d+k(\cos(2\theta)-1)\cos(\theta)\cos(m\phi)$ in polar
coordinate with $(d,k,m)=(6,0.206,2)$. The parameters for the Brusselator
model were chosen to be the same as those used in Fig. 2, indicating that the
pattern is static on a flat plane. We checked the initial condition
sensitivity of the dynamics by comparing two time series of numerical
simulations with slightly different initial conditions; the almost identical
time series in the early stage eventually become distinguished in the late
stage (Fig. 13(b); see details [65]). These results show that chaotic pattern
dynamics can arise due to the surface geometry.

##  VII Discussion

In our previous study, we reported that a Turing pattern that is static on a
flat plane can propagate on curved surfaces [57]. Using the relationship in
Eq. (7), in which the pattern is assumed to move at a constant speed, we
clarified that losses of symmetry in the surface and pattern, such as
reflection asymmetry along the $x$-axis of an axisymmetric surface, are
required for pattern propagation.

In this study, we performed weakly nonlinear analysis as a complementary
approach. This analysis allowed us to explore the generic dynamics of the
Turing pattern on curved surfaces near the bifurcation point. The analysis was
also advantageous in that the surface and pattern symmetries, which were
treated separately in the previous study [57], could be integrated using the
eigenfunctions of the Laplace–Beltrami operator, which are determined by the
surface shape. Because the mode decomposition of the helical and parallel
patterns revealed that static and propagating patterns are described by the
two paired modes, we derived the amplitude equations of two complex variables
corresponding to the two paired modes. The obtained amplitude equations
enabled the classification of pattern dynamics in the vicinity of the Turing
bifurcation depending on the symmetry of the pattern and the surface shape.
For a general axisymmetric surface, the steady-state solution of the amplitude
equations, Eqs. (18) and (19), enables pattern propagation. If a surface has
reflection-symmetry or $2\pi$-periodicity along the $x$-axis, the amplitude
equations depend on the combination of the modes. For helical and parallel
patterns, the amplitude equations drop several terms and are reduced to Eqs.
(30) and (31), for which a propagation solution is no longer possible (see
below for a nontrivial solution). On a flat plane, the amplitude equations are
even simpler (Eqs. (47) and (48)), and the pattern caused by the Turing
instability does not move. Overall, the interaction between the modes is more
intense on general surfaces and leads to the onset of moving patterns. The
propagation conditions given by the relationship between the surface and
pattern symmetry are consistent with previous results.

The present study revealed several points beyond previous results. First, we
found that the periodicity along the $x$-axis is involved in determining the
pattern dynamics. The parallel pattern does not move on the $2\pi$-periodic
surface (Fig. 1), whereas the loss of periodicity causes pattern propagation
(Fig. 8). Second, the pattern can propagate even on a surface with reflection-
symmetry and $2\pi$-periodicity, corresponding to the nontrivial steady-state
solution of the amplitude equations. For such a solution, the pattern profile
is not symmetric. This type of solution is consistent with the propagating
pattern on a sphere determined numerically in the previous study [57]. Third,
a weakly nonlinear analysis revealed that patterns on axisymmetric surfaces
exhibited more complex dynamics than propagation, including oscillating and
chaotic chemical waves. Some of these dynamics were confirmed through
numerical simulations of the RDE for the Brusselator model (Fig. 11). In
particular, motivated by the chaotic solutions in the amplitude equations, we
found a chaotic Turing pattern by direct numerical simulation of Brusselator
model on a deformed spherical surface (see the joint paper [65]). These
findings provide a new perspective on the intricate interplay between surface
characteristics and pattern dynamics.

The present study should be extended to pattern dynamics on general surfaces
without axisymmetry, for which weakly nonlinear analysis is applicable. By
performing numerical simulations of the RDE on several surfaces, we observed
chaotic patterns (Fig. 13) and moving patterns with varying speeds (see Movie
2 [68]). Here, we briefly consider these patterns in terms of the amplitude
equations: without axisymmetry (i.e., no translational invariance along the
$\theta$-axis), amplitude equations can have quadratic terms representing
additional interactions among modes, and more complex dynamics could arise. A
detailed analysis is required in the future.

Pattern dynamics driven by surface curvature are applicable to natural and
engineering systems. By controlling the parameters of the surfaces and
reactions, the pattern can be switched between oscillating, static, and
propagating states. In biological systems, such dynamics can be used for
information transduction, guidance of cellular migration, and positioning of
molecular localization. Analyses using amplitude equations are useful for
predicting pattern dynamics and can also provide a basis for the engineering
of patterns on curved surfaces. In the future, analysis of the stability and
bifurcation of the amplitude equations will further improve our understanding
of pattern dynamics on general curved surfaces, and will also be useful for
dynamics in network systems [72, 73], and on deformable surfaces [74, 75, 76].

Acknowledgements
This study was supported by JST SPRING JPMJSP2108 (to R.N.), JSPS KAKENHI
JP23KJ0641 (to R.N.), JPJSJPR 20191501 (to S.I.), and JST CREST JPMJCR1923,
Japan (to S.I.).

Author contributions
R.N. and S.I. both proposed the research direction, contributed to the
theoretical analysis, and wrote the manuscript. R.N. performed all numerical
simulations.

## References

  * [1] Y. Fuseya, H. Katsuno, K. Behnia and A. Kapitulnik, Nanoscale Turing patterns in a bismuth monolayer, Nature Physics $\bm{17}$, 1031-1036 (2021).
  * [2] T. Aasba, L. Peng, T. Ono, S. Akutagawa, I. Tanaka, H. Murayama, S. Suetsugu, A. Razpopov, Y. Kasahara, T. Terashima et al., Growth of self-integrated atomic quantum wires and junctions of a Mott semiconductor, Science Advances $\bm{9}$, 173-179 (2023).
  * [3] J. Ross, S. C. Müller, and C. Vidal, Chemical waves, Science $\bm{240}$, 460 (1988).
  * [4] J. D. Murray, Mathematical Biology I: An Introduction Vol. 3., (Springer, New York, 2002), J. D. Murray, Mathematical biology II: Spatial models and biomedical applications. Vol. 3., (Springer, New York, 2003).
  * [5] L. Manukyan, S. A. Montandon, A. Fofonjka, S. Smirnov and M. C. Milinkovitch, A living mesoscopic cellular automaton made of skin scales, Nature $\bm{544}$, 173-179 (2017).
  * [6] C. E. Tarnita, J. A. Bonachela, E. Sheffer, J. A. Guyton, T. C. Coverdale, R. A. Long and R. M. Pringle, A theoretical foundation for multi-scale regular vegetation patterns, Nature $\bm{541}$, 398-401 (2017).
  * [7] R. Goyal, M. Jucker, A. Sen Gupta, H. H. Hendon, and M. H. England, Zonal wave 3 pattern in the southern hemisphere generated by tropical convection, Nat. Geosci. $\bm{14}$, 732 (2021).
  * [8] A. Radja, E. M. Horsley, M. O. Lavrentovich, and A. M. Sweeney, Pollen cell wall patterns form from modulated phases, Cell $\bm{176}$, 856 (2019).
  * [9] S. Alonso, M. Bär, and B. Echebarria, Nonlinear physics of electrical wave propagation in the heart: a review, Rep. Prog. Phys. $\bm{79}$, 096601 (2016).
  * [10] J. Bischof, C. A. Brand, K. Somogyi, I. Májer, S. Thome, M. Mori, U. S. Schwarz, and P. Lénárt, A cdk1 gradient guides surface contraction waves in oocytes, Nat. Commun. $\bm{8}$, 1-10 (2017).
  * [11] L. Wettmann and K. Kruse, The Min-protein oscillations in Escherichia Coli: an example of self-organized cellular protein waves, Philos. Trans. R. Soc. Lond. B Biol. Sci. $\bm{373}$, (2018).
  * [12] T. H. Tan, J. Liu, P. W. Miller, M. Tekant, J. Dunkel, and N. Fakhri, Topological turbulence in the membrane of a living cell, Nat. Phys. $\bm{16}$, 657-662 (2020).
  * [13] N. Saito, and S. Sawai, Three-dimensional morphodynamic simulations of macropinocytic cups, iScience $\bm{24}$, 103087 (2021).
  * [14] M. K. McGuire, C. A. Fuller, J. F. Lindner, and N. Manz, Geographic tongue as a reaction-diffusion system, Chaos $\bm{31}$, 033118 (2021).
  * [15] Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis, Nat. Phys. $\bm{17}$, 251-259 (2021).
  * [16] D. Baptista, L. Teixeira, C. van Blitterswijk, S. Giselbrecht, and R. Truckenmüller, Overlooked? underestimated? effects of substrate curvature on cell behavior, Trends Biotechnol. $\bm{37}$, 838 (2019).
  * [17] S. Ehrig, B. Schamberger, C. M. Bidan, A. West, C. Jacobi, K. Lam, P. Kollmannsberger, A. Petersen, P. Tomancak, K. Kommareddy et al., Surface tension determines tissue shape and growth kinetics, Sci. Adv., $\bm{5}$, eaav9394 (2019).
  * [18] B. Schamberger, R. Ziege, K. Anselme, M. B. Amar, M. Bykowski, A. P. G. Castro, A. Cipitria, R. A. Coles, R. Dimova, M. Eder et al., Curvature in biological systems: Its quantification, emergence, and implications across the scales, Adv. Mater., $\bm{35}$, e2206110 (2023).
  * [19] B. Antony, Mechanisms of membrane curvature sensing, Annu. Rev. Biochem. $\bm{80}$, 101-123 (2011).
  * [20] M. Simunovic, G. A. Voth, A. Callan-Jones, and P. Bassereau, When physics takes over: BAR proteins and membrane curvature, Trends Cell Biol. $\bm{25}$, 780 (2015).
  * [21] S. Yang, X. Miao, S. Arnold, B. Li, A. T. Ly, H. Wang, M. Wang, X. Guo, M. M. Pathak, W. Zhao et al., Membrane curvature governs the distribution of Piezo1 in live cells, Nat. Commun. $\bm{13}$, 7467 (2022).
  * [22] L. Pieuchot, J. Marteau, A. Guignandon, T. Dos Santos, I. Brigaud, P. F. Chauvy, T. Cloatre, A. Ponche, T. Petithory, P. Rougerie et al., Curvotaxis directs cell migration through cell-scale curvature landscapes, Nat. Commun. $\bm{9}$, 1-13 (2018).
  * [23] M. Luciano, S.-L. Xue, W. H. De Vos, L. Redondo-Morata, M. Surin, F. Lafont, E. Hannezo, and S. Gabriele, Cell monolayers sense curvature by exploiting active mechanics and nuclear mechanoadaptation, Nature Physics, $\bm{17}$, 1382-1390 (2021).
  * [24] K. Horibe, K. Hironaka, K. Matsushita, and K. Fujimoto, Curved surface geometry-induced topological change of an excitable planar wavefront, Chaos $\bm{29}$, 093120 (2019).
  * [25] V. A. Davydov, V. G. Morozov, and N. V. Davydov, Ring-shaped autowaves on curved surfaces, Phys. Lett. A $\bm{267}$, 326-330 (2000).
  * [26] J. Gomatam and F. Amdjadi, Reaction-Diffusion Equations on a Sphere: Meandering of Spiral Waves, Phys. Rev. E $\bm{56}$, 3913 (1997).
  * [27] H. Yagisita, M. Mimura, and M. Yamada, Spiral wave behaviors in an excitable reaction-diffusion system on a sphere, Physica D $\bm{124}$, 126 (1998).
  * [28] V. A. Davydov, V. G. Morozov, and N. V. Davydov, Critical properties of autowaves propagating on deformed cylindrical surfaces, Phys. Lett. A $\bm{307}$, 265 (2003).
  * [29] F. Kneer, E. Schöll, and M. A. Dahlem, Nucleation of reaction-diffusion waves on curved surfaces, New J. Phys. $\bm{16}$, 053010 (2014).
  * [30] S. Białecki, P. Nałecz-Jawecki, B. Kaźmierczak, and T. Lipniacki, Traveling and standing fronts on curved surfaces, Physica D $\bm{401}$, 132215 (2020).
  * [31] L. R. Gómez, N. A. García, V. Vitelli, J. Lorenzana, and D. A. Vega, Phase nucleation in curved space, Nat. Commun. $\bm{6}$, 6856 (2015).
  * [32] A. R. Singh, T. Leadbetter, and B. A. Camley, Sensing the shape of a cell with reaction diffusion and energy minimization, Proc. Natl. Acad. Sci. U. S. A. $\bm{119}$, e2121302119 (2022).
  * [33] B. Q. Ai, W. J. Zhu, and J. J. Liao, Collective transport of polar active particles on the surface of a corrugated tube, New J. Phys. $\bm{21}$, 093041 (2019).
  * [34] J. P. Brasselet, J. Seade, and T. Suwa, Vector fields on singular varieties (Springer Science & Business Media, 2009).
  * [35] S. Kralj, R. Rosso, and E. G. Virga, Curvature control of valence on nematic shells, Soft Matter $\bm{7}$, 670-683 (2011).
  * [36] L. N. Carenza, G. Gonnella, D. Marenduzzo, G. Negro, and E. Orlandini, Cholesteric Shells: Two-Dimensional Blue Fog and Finite Quasicrystals, Phys. Rev. Lett. $\bm{128}$, 027801 (2022).
  * [37] S. Shankar, M. J. Bowick, and M. C. Marchetti, Topological sound and flocking on curved surfaces, Phys. Rev. X $\bm{7}$, 031039 (2017).
  * [38] F. C. Keber, E. Loiseau, T. Sanchez, S. J. Decamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, Topology and dynamics of active nematic vesicles, Science $\bm{345}$, 1135-1139 (2014).
  * [39] A. M. Turing, The chemical basis of morphogenesis, Phil. Trans. R. Soc. Lond. B $\bm{237}$, 37-72 (1952).
  * [40] S. Kondo, and T. Miura, Reaction-diffusion model as a framework for understanding biological pattern formation, Science $\bm{329}$, 1616-1620 (2010).
  * [41] J. Raspopovic, L. Marcon, L. Russo, and J. Sharpe, Digit patterning is controlled by a Bmp-Sox9-Wnt Turing network modulated by morphogen gradients, Science $\bm{345}$, 566-570 (2014).
  * [42] S. Ishihara, and K. Kaneko, Turing pattern with proportion preservation, J. Theor. Biol. $\bm{238}$, 683-693 (2006).
  * [43] R. G. Plaza, F. Sánchez-Garduño, P. Padilla, R. A. Barrio, and P. K. Maini, The effect of growth and curvature on pattern formation, J. Dyn. Differ. Equ. $\bm{16}$, 1093-1121 (2004).
  * [44] V. Giulio, M. Davide, B. G. Andrew, and O. Enzo, Curvature-driven positioning of Turing patterns in phase-separating curved membranes, Soft Matter $\bm{12}$, 3888-3896 (2016).
  * [45] A. L. Krause, E. A. Gaffney, P. K. Maini, and V. Klika, Modern perspectives on near-equilibrium analysis of Turing systems, Phil. Trans. R. Soc. A. $\bm{379}$, 20200268 (2021).
  * [46] L. Charette, Pattern formation on curved surfaces, Ph.D. Thesis, The University of British Columbia, Canada, (2019).
  * [47] C. Varea, J. L. Aragon, and R. A. Barrio, Turing patterns on a sphere, Phys. Rev. E $\bm{60}$, 4588 (1999).
  * [48] P. C. Matthews, Pattern formation on a sphere, Phys. Rev. E $\bm{67}$, 036206 (2003).
  * [49] M. Núñez-López, G. Chacón-Acosta, and J. A. Santiago, Diffusion-driven instability on a curved surface: spherical case revisited, Braz. J. Phys. $\bm{47}$, 231-238 (2017).
  * [50] D. Lacitignola, B. Bozzini, M. Frittelli, and I. Sgura, Turing pattern formation on the sphere for a morphochemical reaction-diffusion model for electrodeposition, Commun. Nonlinear Sci. Numer. Simul. $\bm{48}$, 484-508 (2017).
  * [51] F. Sánchez-Garduno, A. L. Krause, J. A. Castillo, and P. Padilla, Turing-Hopf patterns on growing domains: the torus and the sphere, J. Theor. Biol. $\bm{481}$, 136-150 (2019).
  * [52] S. S. Liaw, C. C. Yang, R. T. Liu, and J. T. Hong Turing model for the patterns of lady beetles, Phys. Rev. E $\bm{64}$, 041909 (2001).
  * [53] W. Nagata, H. R. Z. Zangeneh, and D. M. Holloway, Reaction-diffusion patterns in plant tip morphogenesis: bifurcations on spherical caps, Bull. Math. Biol. $\bm{75}$, 2346-2371 (2013).
  * [54] S. Nampoothiri, and A. Medhi, Role of curvature and domain shape on Turing patterns, arXiv:1705.02119 (2017).
  * [55] S. Nampoothiri, Effect of geometry on the positioning of a single spot in reaction–diffusion systems, arXiv:1909.06528 (2019).
  * [56] J. R. Frank, J. Guven, M. Kardar, and H. Shackleton, Pinning of diffusional patterns by non-uniform curvature, EPL $\bm{127}$, 48001 (2019).
  * [57] R. Nishide, and S. Ishihara, Pattern Propagation Driven by Surface Curvature, Phys. Rev. Lett. $\bm{128}$, 224101 (2022).
  * [58] G. Nicolis, and I. Prigogine, Self–organization in nonequilibrium systems (Wiley, New York, 1977).
  * [59] B. Peña and C. Pérez-García, Stability of Turing Patterns in the Brusselator Model, Phys. Rev. E $\bm{64}$, 056213 (2001).
  * [60] I. Lengyel, and I. R. Epstein, Modeling of Turing structures in the chlorite-iodide-malonic acid-starch reaction system, Science $\bm{251}$, 650-652 (1991).
  * [61] T. Bánsaǵi Jr., and A. F. Taylor, Helical Turing patterns in the Lengyel-Epstein model in thin cylindrical layers, Chaos $\bm{25}$, 064308 (2015).
  * [62] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. $\bm{65}$, 851 (1993).
  * [63] R. Hoyle, Pattern formation: An introduction to methods (Cambridge University Press, 2006).
  * [64] Y. Kuramoto, Chemical oscillations, waves, and turbulence (Springer Science & Business Media, 2012).
  * [65] R. Nishide, and S. Ishihara, Oscillatory and chaotic pattern dynamics driven by surface curvature.
  * [66] A. K. M. Nazimuddin, M. H. Kabir, and M. O. Gani, Oscillatory wave patterns and spiral breakup in the Brusselator model using numerical bifurcation analysis, J. Comput. Sci. $\bm{62}$, 101720 (2022).
  * [67] A. L. Krause, M. A. Ellis, and R. A. Van Gorder, Influence of curvature, growth, and anisotropy on the evolution of Turing patterns on growing manifolds, Bull. Math. Biol. $\bm{81}$, 759-799 (2019).
  * [68] Supplemental Material.
  * [69] G. Xu, Discrete Laplace–Beltrami operator on sphere and optimal spherical triangulations, Int. J. Comput. Geom. Appl. $\bm{16}$, 75-93 (2006).
  * [70] J. Porter and E. Knobloch, Complex dynamics in the 1:3 spatial resonance, Physica D $\bm{143}$, 138 (2000).
  * [71] C. Kittel and P. McEuen, Introduction to solid state physics (John Wiley & Sons, 2018).
  * [72] H. Nakao, and A. S. Mikhailov, Turing patterns in network-organized activator-inhibitor systems, Nat. Phys. $\bm{6}$, 544-550 (2010).
  * [73] J. van der Kolk, G. García-Pérez, N. E. Kouvaris, M. Á. Serrano, and M. Boguñá, Emergence of Geometric Turing Patterns in Complex Networks, Phys. Rev. X $\bm{13}$, 021038 (2023).
  * [74] G. Toole and M. K. Hurdal, Turing models of cortical folding on exponentially and logistically growing domains, Comput. Math. Appl. $\bm{66}$, 1627 (2013).
  * [75] P. W. Miller, N. Stoop, and J. Dunkel, Geometry of wave propagation on active deformable surfaces, Phys. Rev. Lett. $\bm{120}$, 268001 (2018).
  * [76] N. Tamemoto, and H. Noguchi, Pattern formation in reaction-diffusion system on membrane with mechanochemical feedback, Sci. Rep. $\bm{10}$, 19582 (2020).

Figure 1: Relationships among pattern propagation, pattern profile, and
surface symmetry reported in [57]. The helical pattern propagates on a
reflection-asymmetric surface. The surface shapes are given by the radial
function $r(x)=1.7+0.3\cos(x)+0.05\cos(2x)$ (for the reflection-symmetric
surface) and $r(x)=1.7+0.3\cos(x)+0.05\sin(2x)$ (for the reflection-asymmetric
surface with $-2\pi\leq x<2\pi$). The color map represents $u(x,\theta)$, with
a more intense yellow color indicating a higher chemical concentration.
Figure 2: Propagating pattern on axisymmetric surface. (a) Axisymmetric
surface; radial function $r(x)$ is shown in the right. (b) Dispersion relation
of Brusselator model with parameter set $(a,b,D_{u},D_{v})=(2.0,4.5,0.5,1.8)$.
(c) Static pattern on a flat surface for Brusselator model. (d) Propagating
pattern (left panel) and kymograph along the $\theta$ axis on the pale red
line (right panel). The surface radius is $r(x)=d+k_{1}\cos(x)+k_{2}\sin(2x)$,
where $(d,k_{1},k_{2})=(1.7,0.3,0.05)$. Figure 3: Mode decomposition of the
helical propagating pattern. (a,b) Helical propagating pattern and projected
image on $x$-$\theta$ planes. The parameters are the same as those in Sec.
II.3. (c) Approximated helical propagating pattern constructed using the two
paired modes. (d) Dispersion relation of the Brusselator model. The blue
($\lambda_{1}=2.58$) and orange ($\lambda_{2}=2.32$) points indicate the two
paired modes of the approximation. (e) Eigenfunctions of the two paired modes.
$X_{1}(x)\cos(2\theta),X_{1}(x)\sin(2\theta),X_{2}(x)\cos(2\theta)$, and
$X_{2}(x)\sin(2\theta)$. (f) Time evolution of the amplitudes of the two
paired modes. Real and imaginary parts of $C_{1}$ and $C_{2}$ are shown. (g)
Time evolution of phases $\Theta_{1}$ (blue) and $\Theta_{2}$ (orange) for the
two paired modes.  Figure 4: Pattern consisting of single paired mode. (a,b)
Static pattern and projected image on $x$-$\theta$ plane. We used parameter
set $(a,b,D_{u},D_{v})=(1.7,3.8,0.5,1.8)$ and surface
$r(x)=1.7+0.3\cos(x)+0.05\sin(2x)$. (c) Approximated pattern constructed using
a single paired mode. The Laplace–Beltrami eigenvalue of the mode is
$\lambda=1.92$. (d) The amplitude $C=|C|e^{i\theta}$. (e) Laplace–Beltrami
eigenfunctions $X_{1}\cos(2\theta)$ (left) and $X_{2}\sin(2\theta)$ (right).
Figure 5: Patterns corresponding to the steady-state solutions of Eqs.
(32)–(34). (a–c) Patterns of the Brusselator model on axisymmetric surfaces,
their projections onto the $x$-$\theta$ planes, and the approximated patterns
constructed by eigenmodes. The amplitudes $C_{1}$ (blue) and $C_{2}$ (orange),
and the corresponding eigenfunctions $X_{1}(x)$ and $X_{2}(x)$ are also shown.
The phases of $C_{2}$ are set as $\Theta_{2}=0$. (c) Time evolution of the
pattern and kymograph along the $\theta$ axis on the pale red line. Surface
shapes $r(x)$ are described in the main text. (a) Static pattern on
$2\pi$-periodic surface. $(\lambda_{1},\lambda_{2})=(1.78,1.42)$,
$(a,b,D_{u},D_{v})=(1.77,3.25,0.5,2.5)$. (b) Static pattern on reflection-
symmetric surface. $(\lambda_{1},\lambda_{2})=(1.84,1.86)$,
$(a,b,D_{u},D_{v})=(2.07,3.73,0.5,2.5)$. (c) Propagating pattern on
$2\pi$-periodic surface. $(\lambda_{1},\lambda_{2})=(2.04,1.51)$,
$(a,b,D_{u},D_{v})=(1.62,3.55,0.5,1.7)$. Figure 6:  Patterns corresponding to
the steady-state solutions of the amplitude equations in Eqs. (21)–(23). The
same configuration as Fig. 5. Surface shapes $r(x)$ are described in the main
text. (a) Static pattern on reflection-symmetric surface.
$(\lambda_{1},\lambda_{2})=(1.44,1.97)$,
$(a,b,D_{u},D_{v})=(1.64,3.43,0.5,1.9)$. (b) Static pattern on $2\pi$-periodic
surface. $(\lambda_{1},\lambda_{2})=(1.60,1.36)$,
$(a,b,D_{u},D_{v})=(1.65,3.05,0.5,2.5)$. (c) Propagating pattern on general
surface. $(\lambda_{1},\lambda_{2})=(1.83,1.91)$,
$(a,b,D_{u},D_{v})=(2.09,3.77,0.5,2.5)$. Figure 7: Eigenfunctions and phase
differences of parallel and helical patterns. (a) Parallel static pattern on
$2\pi$-periodic and reflection-symmetric surface. (b) Parallel static pattern
on $2\pi$-periodic surface. (c) Helical static pattern on $2\pi$-periodic and
reflection-symmetric surface. (d) Helical propagating pattern on
$2\pi$-periodic surface. Surfaces are chosen as
$r(x)=1.7+0.3\cos(x)+0.05\cos(2x)$ (a,c) and
$r(x)=1.7+0.3\cos(x)+0.05\sin(2x)$ (b,d). Parameters are set as
$(a,b,D_{u},D_{v})=(2.0,4.5,0.5,1.8)$ (a,b) and $(2.8,5,0.4,2.4)$ (c,d). See
also Supplemental Figures S2–S4 [68].  Figure 8: Propagation of a parallel
pattern on a $2\pi$-aperiodic surface. (a) Surface shape
$r(x)=1.7+0.3\cos(x)+0.05\sin(2x)+0.1\cos(x/2)$. (b) Pattern projected onto
the $x$-$\theta$ plane. (c) Time evolution of patterns on the surface and
kymograph along the $\theta$-axis, represented by a pale red line. The
parameters are set as $(a,b,D_{u},D_{v})=(2.8,5,0.4,2.4)$.  Figure 9:
Comparison of steady solutions for the Brusselator model and amplitude
equations. (a)–(d) Trajectories $(\eta,\xi,\alpha)$ for the amplitude
equations (blue) and the Brusselator model (orange) are shown in the left
panels, with steady-state solutions indicated by large points. Patterns on
surfaces and their projected images onto $x$-$\theta$ planes are also shown
for the amplitude equations and the original Brusselator model. Parameters are
set as (a) $r(x)=1.7+0.18\cos(x)+0.03\cos(2x)$,
$(a,b_{0},D_{u},D_{v})=(2.29,4.86,0.5,1.8)$ and $\nu=0.05$. (b)
$r(x)=1.95+0.21\cos(x)-0.06\cos(2x)+0.015\cos(x/2)$,
$(a,b_{0},D_{u},D_{v})=(1.64,3.41,0.5,1.9)$ and $\nu=0.02$. (c)
$r(x)=1.7+0.18\cos(x)+0.03\sin(2x)$,
$(a,b_{0},D_{u},D_{v})=(2.29,4.86,0.5,1.8)$ and $\nu=0.05$. (d)
$r(x)=1.9+0.21\cos(x)+0.06\sin(2x)$,
$(a,b_{0},D_{u},D_{v})=(1.62,3.54,0.5,1.8)$ and $\nu=0.01$.  Figure 10:
Approximation accuracy of amplitude equations for Brusselator model. (a)
Amplitudes $\eta$ and (b) $\xi$, (c) phase difference $\alpha$, and (d)
velocity $\omega$ are plotted against bifurcation parameter $\nu$. The blue
and orange lines represent the amplitude equations and the original
Brusselator model, respectively. Figure 11: Limit cycle solution. (a) Surface
$r(x)=1+0.15\sin(2x)+0.4\sin(3x)$ (b) Dispersion relation and Laplace–Beltrami
eigenvalues corresponding to the two paired modes (orange and green points).
(c) Trajectory of the limit cycle solution $(\eta,\xi,\alpha)$ (left) and time
series of $\eta,\xi$ and $\alpha$ (right). The blue and orange lines represent
the amplitude equations and the original Brusselator model, respectively. (d)
Snapshots of the limit cycle solutions and kymographs along the $\theta$-axis
represented by pale red lines for the amplitude equations (left) and
Brusselator model (right). See also Movie 1 [68]. Figure 12: Chaotic dynamics
in amplitude equations. (a) Trajectory of chaotic solution at $b_{3}=-0.09$.
(b) Bifurcation diagram (top) and Lyapunov exponent (bottom) against $b_{3}$.
For the bifurcation diagram, the value of $\eta$ on a Poincaré section is
determined by $d\eta/dt=0$ was plotted for each value of $b_{3}$. Figure 13:
Chaotic dynamics of Turing pattern on a deformed sphere. (a) Time evolution of
Turing pattern. (b) Time series of chemical concentrations $U=u-a$ at a point
of the surface. The orange and light blue lines represent two different time
series with slightly different initial conditions. See the joint paper [65]
for details.