# Oscillating Turing patterns, chaos and strange attractors in a reaction-diffusion system augmented with self-and cross-diffusion terms

Benjamin Aymard

► To cite this version:

Benjamin Aymard. Oscillating Turing patterns, chaos and strange attractors in a reaction-diffusion system augmented with self-and cross-diffusion terms. 2025. hal-04887186

HAL Id: hal-04887186

https://inria.hal.science/hal-04887186v1

Preprint submitted on 14 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Oscillating Turing patterns, chaos and strange attractors in a reaction-diffusion system augmented with self- and cross-diffusion terms

Benjamin Aymard

MathNeuro Team,
Inria branch at the University of Montpellier,
Montpellier, France

## Abstract

In this article we introduce an original model in order to study the emergence of chaos in a reaction diffusion system in the presence of self- and cross-diffusion terms. A Fourier Spectral Method is derived to approximate equilibria and orbits of the latter. Special attention is paid to accuracy, a necessary condition when one wants to catch periodic orbits and to perform their linear stability analysis via Floquet multipliers. Bifurcations with respect to a single control parameter are studied in four different regimes of diffusion: linear diffusion, self-diffusion for each of the two species, and cross-diffusion. Key observations are made: development of original Turing patterns, Hopf bifurcations leading to oscillating patterns and period doubling cascades leading to chaos. Eventually, original strange attractors are reported in phase space.

Keywords: Self-diffusion; Cross-diffusion; Turing instability; Hopf bifurcation; Floquet theory; Period doubling cascade; Chaos; Strange attractors.

## 1. Introduction

When does a deterministic dynamical system become unpredictable? In his 1814 memoir on probabilities [1], Pierre-Simon de Laplace hypothetised that one could predict the future of a dynamical system given the exact knowledge of positions and impulsion of its components at an initial time. Later, James
Clerk Maxwell observed in 1876 [2] that for a class of systems, an error on the initial condition could dramatically modify its outcome. Eventually, Henri Poincaré's work in celestial mechanics from 1892 [4] led him to discover the unpredictability of trajectories in the three body problem. In doing so, he provided fundamental tools, of what would later become the chaos theory. Beyond its philosophical aspect, and its origin in celestial mechanics, chaos theory finds practical applications in everyday life, in problems ranging from meteorology [10] to electrical engineering, passing by chemical reactions [12], plasma physics [9], laser physics [23] and mechanics [22]. Understanding chaotic systems is a challenging task, with numerous potential applications.

Using the tools of bifurcation theory [28, 30], several mechanisms explaining the development of chaos have been identified. The key idea of this approach is to evaluate the sequence of bifurcations (loss of stability of equilibria or periodic orbits) while a control parameter is varied. Among the identified scenarios is the period doubling cascade. Discovered in 1976 in numerical experiments on the logistic map (see [11]), then observed in numerous experiments, such as in electronic circuits in 1981 by Linsay [18], or in mercury convection rolls in 1982 by Libchaber [19], it consists of a sequence of period doubling of periodic orbits. After each period doubling, the distance between consecutive bifurcation diminishes, and eventually, the period tends to infinity for a finite value of the control parameter: the orbit becomes aperiodic. Careful measurements of distances between consecutive bifurcation values led to the discovery of the Feigenbaum constant [15, 16], revealing a concept of universality of the period doubling cascade. It was then observed in many other models, in particular in the Rössler model for chemical reaction [12], the Mackey-Glass system in blood cells population dynamics [13], neuron models of Fitzhugh-Nagumo type [20], the Kuramoto-Sivashinsky system modeling flame front propagation [21], the Chua circuit [24], the forced Dufing oscillator [5], among others. Interested reader may find a historical note about period doubling in [34].

However, even if many roads to chaos have been identified, the question of understanding how does the presence of self- and/or cross-diffusion affect the
development of chaos in reaction diffusion systems, remains unclear. This work, which aims precisely at answering this question, stands at the crossroads of three domains that appeared during the 20th century. First, the domain of morphogenesis, initiated by the seminal paper of Alan Turing [8], studying the appearance of periodically spatial patterns, as a result of an instability created by the interaction between linear diffusion, and nonlinear reaction. Second, the domain of spatio-temporal chaos, notably following the works of Kuramoto on diffusion induced chaos, where, for instance, the interaction between a Hopf bifurcation and a Turing instability can create oscillating patterns, possibly becoming chaotic [7, 14]. Third, the domain of self- and cross-diffusion, motivated by questions in population dynamics and in chemistry, notably the works of Shigesada Kawasaki Terramoto [17] and Vanag and Epstein [27], and more recently in [35, 39].

There is therefore a need for a simple model, with bifurcations driven by a single parameter, gathering those three elements, for which we could study in a continuous way the development of chaos. Due to the presence of nonlinearities on both reaction and diffusion terms, this kind of system is, in general, not analytically solvable, except in special cases. Therefore, reliable numerical methods [29] have to be designed in order to study them numerically. To this end, high accuracy is a necessary condition when one wants to catch periodic orbits, and to perform their linear stability analysis via Floquet multipliers [3].

The outline of the paper is the following. First, an original model is introduced, containing self- and cross-diffusion, with bifurcations controled by a single parameter. Together with the model, a Fourier spectral method is derived, and theoretical convergence guarantees are provided. Second, bifurcations, of both equilibria and periodic orbits, are studied in four regimes: linear diffusion, activator self-diffusion, inhibitor self-diffusion, and cross-diffusion. Finally, the road to chaos is identified in each case, and the strange attractors that appear in the chaotic regime, are described.
## 2. Model and method

### 2.1. Model

In [26, 33] the authors introduce a model containing a richness of behavior, while being relatively simple, as the bifurcations are driven by one control parameter C, making it, as they mention, a good laboratory for studying oscillation and chaos in two species reaction diffusion models. We will refer to this model as the classical BVAM model. Along those lines, our aim is to introduce a generalization of their model, by introducing non linear diffusion terms. Let us consider the dynamics of two species, $u_1$ and $u_2$, on a domain $\Omega = [-L_x, L_x]$, described by the reaction diffusion system:

$$ 
\begin{aligned}
\frac{\partial u_1}{\partial t} &= \Delta(d_1 u_1 + d_{11} u_1^3 + d_{12} u_2^2 u_1) + \eta(u_1 + au_2 - Cu_1 u_2 - u_1 u_2^2), \\
\frac{\partial u_2}{\partial t} &= \Delta(d_2 u_2 + d_{22} u_2^3 + d_{12} u_1^2 u_2) + \eta(bu_2 + Hu_1 + Cu_1 u_2 + u_1 u_2^2),
\end{aligned}
\quad (1) $$

with $\eta$, $a$, $b$, $C$ and $H$ reaction parameters, and $d_i, d_{ii}, d_{ij}$ for $i, j = 1, 2, i \neq j$ diffusion parameters. Let us remark that, when $d_{11} = d_{22} = d_{12} = 0$, the model (1) reduces to the original BVAM model. The problem is mathematically closed by adding initial conditions and periodic boundary conditions on the domain $\Omega$. In their original work [33], the authors considered Dirichlet boundary conditions. Our choice of periodic boundary conditions here is motivated by two reasons: this reduces the importance of the central point of the domain, later used to define a phase space, and it allows to easily use Fourier analysis [25]. From [38] we know that models of the form (1) follow an energy law of the form:

$$ \frac{d}{dt} E(t) = - \sum_{i=1}^{2} \|\nabla \mu_i\|_{L^2(\Omega)}^2 + \sum_{i=1}^{2} (R_i, \mu_i), \quad (2) $$

with $E$ defined by:

$$ E(t) = \int_{\Omega} \left( d_1 \frac{u_1^2}{2} + d_2 \frac{u_2^2}{2} + d_{11} \frac{u_1^4}{4} + d_{22} \frac{u_2^4}{4} + d_{12} \frac{u_1^2 u_2^2}{2} \right) dx, \quad (3) $$

with reaction terms $R_1, R_2$ defined by:

$$ 
\begin{aligned}
R_1 &= \eta(u_1 + au_2 - Cu_1u_2 - u_1u_2^2), \\
R_2 &= \eta(bu_2 + Hu_1 + Cu_1u_2 + u_1u_2^2),
\end{aligned}
\quad (4) $$
and chemical potentials $\mu_1, \mu_2$ defined by:

$$ 
\begin{aligned}
\mu_1 &= d_1 u_1 + d_{11} u_1^3 + d_{12} u_2^2 u_1, \\
\mu_2 &= d_2 u_2 + d_{22} u_2^3 + d_{12} u_1^2 u_2.
\end{aligned} 
\qquad (5) 
$$

Due to the complexity of the model, analytical solutions may only exist in special cases; therefore, an efficient numerical method is needed in order to approximate orbits of this dynamical system.

## 2.2. Fourier Spectral Method

### 2.2.1. Spatial approximation

Using periodicity assumption, one may consider the Fourier series (in space):

$$ u_i(x,t) = \sum_{k=-\infty}^{\infty} (\hat{u}_i)_k(t) e^{ikx/|\Omega|}, \quad i = 1, 2 \qquad (6) $$

$$ \mu_i(x,t) = \sum_{k=-\infty}^{\infty} (\hat{\mu}_i)_k(t) e^{ikx/|\Omega|}, \quad i = 1, 2 \qquad (7) $$

$$ R_i(x,t) = \sum_{k=-\infty}^{\infty} (\hat{R}_i)_k(t) e^{ikx/|\Omega|}, \quad i = 1, 2, \qquad (8) $$

with Fourier coefficients defined by:

$$ (\hat{u}_i)_k(t) = \frac{1}{|\Omega|} \int_{\Omega} u(x,t) e^{-ikx/|\Omega|} dx, \quad i = 1, 2 \qquad (9) $$

$$ (\hat{\mu}_i)_k(t) = \frac{1}{|\Omega|} \int_{\Omega} \mu(x,t) e^{-ikx/|\Omega|} dx, \quad i = 1, 2 \qquad (10) $$

$$ (\hat{R}_i)_k(t) = \frac{1}{|\Omega|} \int_{\Omega} R(x,t) e^{-ikx/|\Omega|} dx, \quad i = 1, 2. \qquad (11) $$

Injecting those expressions in the model (1), and using (4), (5), one may write:

$$ \frac{\partial \hat{u}_i}{\partial t} = -k^2 \hat{\mu}_i + \hat{R}_i, \quad i = 1, 2. \qquad (12) $$

In Fourier space, the set of partial differential equations has become a set of ordinary differential equations, easier to solve. However, except in simple cases, the Fourier series have to be evaluated numerically. To this end, we use the Fast Fourier Transform algorithm [25]. Given a discretization of the domain $\Omega = [-L_x, L_x]$ by a step:

$$ \Delta x = \frac{2L_x}{N}, $$
leading to collocation points:

$$x_j = j\Delta x \text{ with } j = -\frac{N}{2}, \dots, \frac{N-1}{2}.$$

we compute an approximation of (11), called the discrete Fourier transform:

$$\tilde{\mu}_k = \frac{1}{N} \sum_{j=0}^{N-1} \mu_j e^{-ikx_j},$$

Frequencies $f_m$ are evaluated as:

$$f_m = m\Delta f$$

with:

$$\Delta f = \frac{1}{N\Delta x} = \frac{1}{2L_x}.$$

Eventually, wave numbers $k$ and frequencies $f$ are related as:

$$k = 2\pi f.$$

**Remark 1.** Several conventions coexist for the FFT algorithm. The results are independent from the used convention, however, it is crucial to respect the FFT convention of the used library when one defines $k$, in order to be consistent during all the FFT and inverse FFT operations.

**Remark 2.** The reader may note that nonlinear terms are first evaluated in real space, then transformed in Fourier space.

The spatial approximation verifies an important convergence property, that we will now demonstrate. Let us first recall the definition of a type of Sobolev space [25]:

$$H^m = \{f \in L^2, (1 + |k|^2)^{m/2} \hat{f} \in L^2\}$$

**Proposition 1.** Let us assume that model (1) admits equilibrium solutions:

$$-\Delta\mu_i(u_1, u_2) = R_i(u_1, u_2). \qquad (13)$$

Let us further assume that those solutions are bounded:

$$\|u_i\|_{L^2} \le C.$$

and smooth enough such that $u_i \in H^m$. Then, for $N$ large enough, the following inequality holds:

$$\|\mu_i - P_N(\mu_i)\|_{L^2}^2 \le \frac{\|R_i\|_{H^m}^2}{(1 + N^2)^m N^4}. \qquad (14)$$
**Remark 3.** The assumptions result from observations of numerous simulations with various parameters. However, the proof of existence, regularity and boundedness of solutions remain open.

*Proof.* By application of Fourier synthesis theorem [25], one gets:

$$ 
\begin{aligned}
\mu_i(x) &= \sum_{k=-\infty}^{\infty} \hat{\mu}_{ik} e^{ikx} \\
&= \underbrace{\sum_{k=-N/2}^{N/2} (\tilde{\mu}_i)_k e^{ikx}}_{P_N(\mu)} + \underbrace{\sum_{k=-N/2}^{N/2} ((\hat{\mu}_i)_k - (\tilde{\mu}_i)_k) e^{ikx}}_{A_N(\mu)} + \underbrace{\sum_{|k|>N/2} (\hat{\mu}_i)_k e^{ikx}}_{T_N(\mu)}
\end{aligned}
$$

with $P_N$ the interpolation polynomial, $A_N$ the aliasing error, and $T_N$ the truncation error. We then get, using orthogonality of the Fourier basis:

$$ \|\mu_i - P_N(\mu_i)\|_{L^2}^2 = \|A_N(\mu_i)\|_{L^2}^2 + \|T_N(\mu_i)\|_{L^2}^2. $$

For the aliasing error $A_N$, let us recall the classical argument, relying on periodicity of the complex exponential, and the Fourier series evaluated at collocation points:

$$ \mu_i(x_j) = \sum_{k=-\infty}^{\infty} (\hat{\mu}_i)_k e^{ikx} = \sum_{k=0}^{N-1} \left( \sum_{q=-\infty}^{\infty} (\hat{\mu}_i)_{(k+qN)} \right) e^{ikx}, $$

then, by identification, one gets:

$$ (\hat{\mu}_i)_k - (\tilde{\mu}_i)_k = \sum_{q=-\infty, q \neq 0}^{\infty} (\hat{\mu}_i)_{(k+qN)}. $$

For $N$ large enough, as $\mu_1, \mu_2$ are regular, the aliasing error is null. For the truncation error, one may expand, using equation (12):

$$ \|T_N(\mu_i)\|_{L^2}^2 = \sum_{|k|>N/2} |(\hat{\mu}_i)_k|^2 = \sum_{|k|>N/2} \frac{|(\hat{R}_i)_k|^2}{k^4}. $$

Then, multiplying numerator and denominator by the same constant, then by using the hypothesis on regularity, and finally the Parseval identity, one gets:

$$ 
\begin{aligned}
\|T_N(\mu_i)\|_{L^2}^2 &= \sum_{|k|>N/2} \frac{(1+k^2)^m |\hat{R}_{ik}|^2}{(1+k^2)^m k^4} \le \frac{1}{(1+N^2)^m N^4} \sum_{|k|>N/2} (1+k^2)^m |\hat{R}_{ik}|^2 \\
&\le \frac{\|R_i\|_{H^m}^2}{(1+N^2)^m N^4}.
\end{aligned}
$$

As $u_1, u_2$ are bounded, by hypothesis, and as $R_1, R_2$ are polynomials in $u_1, u_2$, then for all $m$, $R_i \in H^m$. $\square$
### 2.2.2. Temporal approximation

To approximate the temporal part, the classical Runge-Kutta method of fourth order (RK4) is used ([31]) Defining a time step $\Delta t$, and denoting by $\hat{u}_i^n$ the approximation of $\hat{u}_i$ at time $t^n = n\Delta t$, the RK4 scheme reads:

$$ \hat{u}_i^{n+1} = \hat{u}_i^n + \frac{\Delta t}{6} (\hat{F}(U_1) + 2\hat{F}(U_2) + 2\hat{F}(U_3) + \hat{F}(U_4)), \quad (15) $$

with:

$$ U_1 = (\hat{u}_1, \hat{u}_2), $$

$$ U_2 = U_1 + \frac{\Delta t}{2} \hat{F}(U_1), $$

$$ U_3 = U_1 + \frac{\Delta t}{2} \hat{F}(U_2), $$

$$ U_4 = U_1 + \Delta t \hat{F}(U_3). $$

## 3. Patterns and oscillations

We study the bifurcations of the system with respect to parameter $C$, in the presence of four different kind of diffusion: a classical linear diffusion (seen as a witness test, the original BVAM model), a self-diffusion on each of the two species, and a cross-diffusion between the two species. Parameters are gathered in 1. Let us recall that in their original article, [33] considered $\Delta x = 0.2$ and $\Delta t = 0.01$.

<table><caption>Table 1: Parameters</caption><thead><tr><th>Reaction</th><th>$L_x = 5, H = 3, \eta = 1, a = -1, b = -3/2$</th></tr></thead><tbody><tr><td>Linear diffusion</td><td>$d_1 = 0.08, d_2 = 1, d_{11} = 0, d_{12} = 0, d_{22} = 0$</td></tr><tr><td>Self-diffusion on $u_1$</td><td>$d_1 = 0.08, d_2 = 1, d_{11} = 0.07, d_{12} = 0, d_{22} = 0$</td></tr><tr><td>Self-diffusion on $u_2$</td><td>$d_1 = 0.08, d_2 = 1, d_{11} = 0, d_{12} = 0, d_{22} = 0.05$</td></tr><tr><td>Cross-diffusion</td><td>$d_1 = 0.08, d_2 = 1, d_{11} = 0, d_{12} = 0.02, d_{22} = 0$</td></tr><tr><td>Spectral method (periodic orbits)</td><td>$N = 300, \Delta t = 8e^{-5}$</td></tr><tr><td>Spectral method (strange attractors)</td><td>$N = 500, \Delta t = 2e^{-5}$</td></tr><tr><td>Newton-Krylov method</td><td>$N = 1000, \epsilon = 1e^{-10}$</td></tr></tbody></table>
### 3.1. Bifurcation of equilibria

#### 3.1.1. Equilibria

Equilibria of dynamical system (1), denoted by $(\bar{u}_1, \bar{u}_2)$, are defined as stationary solutions, verifying:

$$ \frac{\partial}{\partial t} \begin{pmatrix} \bar{u}_1 \\ \bar{u}_2 \end{pmatrix} = \begin{pmatrix} \Delta\mu_1(\bar{u}_1, \bar{u}_2) + R_1(\bar{u}_1, \bar{u}_2) \\ \Delta\mu_2(\bar{u}_1, \bar{u}_2) + R_2(\bar{u}_1, \bar{u}_2) \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \quad (16) $$

Homogeneous steady states are spatially constant equilibria. With the choice of parameters given in Table 1, BVAM have showed that only one solution exist, given by $(\bar{u}_1 = \bar{u}_2 = 0)$. In the presence of sufficiently high diffusion, a Turing instability appears, creating non homogeneous states, that can not be analytically computed in general. In order to evaluate equilibria (16), we use the Newton-Krylov method [32] to find the roots of the equivalent system:

$$ \begin{cases} \mathcal{F}^{-1}(\hat{F}_1(u_1, u_2)) = 0, \\ \mathcal{F}^{-1}(\hat{F}_2(u_1, u_2)) = 0. \end{cases} \quad (17) $$

until a given tolerance $\epsilon$ is reached, starting from an initial guess of the form:

$$ u_i(x) = A_i \cos(n_i x), \quad i = 1, 2. $$

with $n_i$ a parameter, tuning the wave number of the initial guess.

#### 3.1.2. Stability of equilibria

Linear stability analysis consists of studying the stability of equilibria, by studying the flow of linear perturbation around steady states. Given a perturbation $(\delta u_1, \delta u_2)$ around an equilibrium $(\bar{u}_1, \bar{u}_2)$, the linearized flow around a steady state reads:

$$ \frac{\partial}{\partial t} \begin{pmatrix} \delta u_1 \\ \delta u_2 \end{pmatrix} \approx L(\bar{u}_1, \bar{u}_2) \begin{pmatrix} \delta u_1 \\ \delta u_2 \end{pmatrix} $$

with:

$$ L(\bar{u}_1, \bar{u}_2) = \begin{pmatrix} \frac{\partial R_1}{\partial u_1}(\bar{u}_1, \bar{u}_2)\delta u_1 + \Delta \left( \frac{\partial \mu_1}{\partial u_1}(\bar{u}_1, \bar{u}_2)\delta u_1 \right) & \frac{\partial R_1}{\partial u_2}(\bar{u}_1, \bar{u}_2)\delta u_2 + \Delta \left( \frac{\partial \mu_1}{\partial u_2}(\bar{u}_1, \bar{u}_2)\delta u_2 \right) \\ \frac{\partial R_2}{\partial u_1}(\bar{u}_1, \bar{u}_2)\delta u_1 + \Delta \left( \frac{\partial \mu_2}{\partial u_1}(\bar{u}_1, \bar{u}_2)\delta u_1 \right) & \frac{\partial R_2}{\partial u_2}(\bar{u}_1, \bar{u}_2)\delta u_2 + \Delta \left( \frac{\partial \mu_2}{\partial u_2}(\bar{u}_1, \bar{u}_2)\delta u_2 \right) \end{pmatrix} \quad (18) $$
Using FFT, one may approximate operator $L$ in finite dimension. Indeed, as the Fourier Transform is linear, one can evaluate it on a basis, to get a matrix form:

$$F_N = \hat{I}_n.$$

Using this form allows to express Laplace operator as:

$$\Delta_N = F_N^{-1} \times K \times F_N.$$

with $K$ a diagonal matrix, such that $K_{mm} = k_m^2$ and $K_{mn} = 0$ for $m \neq n$. Replacing $\Delta$ operator in (18) by $\Delta_N$ defines a matrix, from which stability is evaluated by computing its spectrum.

### 3.1.3. Results

Using a parameter continuation strategy [30] on parameter $C$, equilibria and their stability have been computed in four regimes: linear diffusion, self-diffusion on $u_1$, then on $u_2$, and eventually cross-diffusion. Results are displayed on Figure 1. For $C$ varying from $-0.5$ to $-1.5$ in 100 steps, steady states are computed, then the corresponding energy (3) is evaluated. Evaluating the $L$ matrix (18) on the steady states, stability of equilibria is also evaluated. Let us note here that a threshold of $1e-3$ is considered on the real part of eigenvalue, to avoid numerical false positives. For each of the four cases, the energy curve grows non linearly as parameter $C$ decreases. Equilibrium is stable at first (continuous line), as the eigenvalues of $L$ are all with a strictly negative real part. Then, stability is lost (dotted line) through a Hopf bifurcation [7] (indicated as a red dot), when a pair of complex eigenvalues, with non null imaginary part, crosses the imaginary axis (their real part becomes strictly positive). Even though the behaviour is qualitatively comparable, the presence of self- and/or cross- diffusion dramatically changes the value of the parameter of Hopf bifurcation.

## 3.2. Bifurcation of periodic orbits

### 3.2.1. Periodic orbits

Let us denote by $\phi^t(X_0)$ the orbit, at time $t$, starting from $X_0 = (u_1(t = 0, x), u_2(t = 0, x))$ at time $t = 0$, solution of the system (1). Periodic orbits are
<figure><img src="image_2.png" /><figcaption>(a) Linear diffusion</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Self-diffusion on $u_1$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Self-diffusion on $u_2$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(d) Cross-diffusion</figcaption></figure>

Figure 1: Bifurcation of equilibria. Steady states (16) of model (1), are computed using (17), and the corresponding energy (3) is displayed, with respect to bifurcation parameter $C$. A Hopf bifurcation appears in each case, when a stable equilibrium (continuous line), becomes unstable (dotted line), as a couple of complex eigenvalues with non negative real part appear in the spectrum (18). Past the Hopf bifurcation, the Turing patterns starts to oscillate in time, creating a periodic orbit.
orbits $\bar{X}$ for which there exists a period $T > 0$ such that:

$$ \forall t, \phi^T(\bar{X}(t)) = \bar{X}(t). \qquad (19) $$

Computing a periodic orbit consists of finding a period $T$, and a corresponding limit cycle $\bar{X}$, solving equation (19). However, this system is not closed. First, $T$ is unknown a priori, and, even for a fixed $T$, there exists a continuum of solution, parametrized by a shift of time $\sigma$. In order to close the problem, a phase condition has to be added: this is the basis of the shooting method [30]. The classical method consists of minimizing the distance $\rho$ between the periodic orbit $\bar{X}$ and the reference solution $\tilde{X}$, including a potential time shift $\sigma$:

$$ \rho(\sigma) = \int_{0}^{T} \|\tilde{X}(\tau + \sigma) - \bar{X}(\tau)\|^2 d\tau. $$

Considering a variation of $\sigma$ leads to the first order condition, cancelling the derivative:

$$ \rho'(\sigma) = 2 \int_{0}^{T} (\tilde{X}(\tau + \sigma) - \bar{X}(\tau)).F(\tilde{X}(\tau))d\tau $$

In our case, we consider, around $\tau = 0$, with $\sigma = T$:

$$ \rho'(T) \approx (\phi^T(\bar{X}(0)) - \bar{X}(0)).F(\tilde{X}(0)). $$

The latter approximation is rather coarse. However, it is fast to evaluate, and works well in practice. On our system, at least, it seems to increase the convergence of the solver. Let us assume that we know a reference solution $\tilde{X} = (\tilde{u}_1, \tilde{u}_2)$. In practice, this will be the solution found using the former parameter value in a continuation loop, and the steady state before the Hopf bifurcation for the initial one. The periodic orbit problem may then be solved using a Newton-Krylov method on the system:

$$ \phi^T(\bar{X}(t)) - \bar{X}(t) = 0, \qquad (20) $$

$$ (\phi^T(\bar{X}(0)) - \bar{X}(0)).F(\tilde{X}(0))) = 0. \qquad (21) $$

For our tests, we have set up a tolerance of $5e^{-4}$ on the overall residual.
### 3.2.2. Floquet theory and orbit stability

Once the periodic orbit $\bar{X}$ is known, together with its period $T$, comes the question of its stability. The classical tool to do so is the concept of Poincaré map, sending points from a section in the phase space (called the Poincaré section), to itself, after a period $T$. In particular, periodic orbits are fixed points of the Poincaré map. Monodromy matrices $M$ (from the Greek "monos" and "dromos", "one road" matrix) are linear approximations of the Poincaré map, defined as:

$$M = \frac{\partial \phi^T(\bar{X}(0))}{\partial x}. \qquad (22)$$

In practice, we will evaluate $M = (c_1, \dots, c_{2N})$ by columns, by a finite difference approximation:

$$c_i \approx \frac{\phi^T(\bar{X}(0) + h e_i) - \phi^T(\bar{X}(0))}{h}.$$

In our applications, we set $h = 0.001$. By application of the Floquet theorem [3], we know that at each eigenvalue $\mu$ of the monodromy matrix, called a Floquet multipliers, corresponds a mode $X_\mu$, such that:

$$\dot{X}_\mu = \mu X_\mu. \qquad (23)$$

Therefore, a mode is stable if and only if the absolute value of the Floquet multiplier $\mu$ is strictly smaller than 1. Stability may be lost by a fold bifurcation ($\mu$ exit the unit circle at 1), a period doubling bifurcation ($\mu$ exits the unit circle at -1), leading to the creation of a cycle of double period, or a Neimark Sacker bifurcation ($\mu$ exits the unit circle at $e^{i\theta}$), leading to the creation of a quasiperiodic cycle (2 periods coexist, but are not commensurate).

### 3.2.3. Results

Using a numerical continuation approach [30], we study the orbits in each of the four cases, starting from the Hopf bifurcation point, where a limit cycle is created. For each case, we display the variation of density of the first species $u_1(t)$, the corresponding variation of its Fourier transform in space $\frac{2}{N}|\mathcal{F}(u - \bar{u})|$, the evolution of energy (3), the Fourier spectrum of the signal of the energy
$\frac{2}{N}|\mathcal{F}(E - \bar{E})|$, the Floquet multipliers when the cycle becomes unstable, and
the periods with respect to the bifurcation parameter.

**Remark 4.** With this choice of convention in the plot of the Fourier transform, a signal of the form $E(t) = A \cos(2\pi f t)$ (resp. $u(x) = A \cos(kx)$) will appear as a peak on the Fourier domain, centered on $f = n$ (resp $k = n$), with height $A$, and therefore will be easily interpretable. Removing the mean values allows to define a refined color bar.

On Figure 2, we consider a variation of $C$ with step $\delta C = -0.01$. We observe
that, in each of the four cases, passed the value of the Hopf bifurcation, a periodic
orbit appears, in perfect agreement with the value computed in the previous
section (see Figure 1). Starting from the unstable steady state, and after an
initial simulation (to be sure to achieve convergence), we use this periodic orbit
as initial condition of our continuation strategy. As initial guess, we set $T = 3$,
according to observations. The Turing patterns, once steady, become unstable,
and start to oscillate in time. For each of the four cases, a period doubling
bifurcation is observed, when the value of parameter $C$ is decreased, as a Floquet
multiplier, in red, crosses the unit circle at $-1$. We notice a peak around $k = 2$
for the density in Fourier domain, leading to a spatial frequency of $f = \frac{2\pi}{k} \approx \pi$,
in agreement with the distribution of density in space.

On Figure 3, we start from the unstable periodic orbit, and we simulate the
orbit until the periodic doubling is observed. We use this orbit as an initial
condition for the shooting method with an initial guess of $T$ as twice the value
of the previous limit cycle. We consider a smaller step of $\delta C = -0.001$, as the
previous step was not refined enough to ensure convergence between consecutive
values of bifurcation parameter. We continue the orbits until they become
unstable. Interestingly, the scenario changes according to the case. For linear
diffusion, self diffusion on $u_1$, and cross-diffusion, a second period doubling is
observed. For the self diffusion on $u_2$, a fold bifurcation appears.

**Remark 5.** In order to check the convergence of numerical results, we have tried $N = 300$ and $N = 500$ in the former computations. We have observed, qualitatively, the same results.
<figure><img src="image_6.png" /><figcaption>Figure 2: Bifurcation of periodic orbits, starting from the Hopf bifurcation (Figure 1). Periodic orbits are computed using the shooting method (21). Density and its Fourier transform (in space) are computed with respect to time, and one full cycle is displayed on lines 1 and 2. Energy (3) and its Fourier transform (in time), on one full cycle, are displayed on lines 3 and 4. Period with respect to control parameter C is displayed on line 5. Floquet multipliers (23) are displayed on line 6, at $C = -1.39, -1.43, -1.48, -1.36$. Eventually, all periodic orbits lose their stability via a period doubling bifurcation ($\mu$ exits the unit circle via $-1$).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3: Bifurcation of periodic orbits, starting from the first period doubling bifurcation (Figure 2). Floquet multipliers are displayed on line 6, at $C = -1.471, -1.482, -1.599, -1.486$. Eventually, all periodic lose stability, either by period doubling ($\mu$ exits the unit circle via -1) or by a fold ($\mu$ exits via 1).</figcaption></figure>
# 4. Chaos and strange attractors

## 4.1. Phase space

In the original BVAM article [33], authors consider attractors in phase space of the form:

$$A(t) = \{(u_1(x = 0, t), u_2(x = 0, t)|t \in [\tau, T]\}. \quad (24)$$

In this work, we also monitor a phase space in terms of energy, by defining:

$$A_E(t) = \left\{ E(t), \frac{d}{dt} E(t) | t \in [\tau, T] \right\}. \quad (25)$$

This attractor is more global, in the sense that, contrarily to the previous one, it is not attached to a specific space coordinate ($x = 0$), but only to a time $t$. In addition, we color the curves with respect to a normalized time (sending $[0, T]$ onto $[0, 1]$); a blue color indicates a small time, while a yellow one a large time.

## 4.2. Road to chaos

On Figure 3, we observed that, for each of the four regimes of diffusion, the periodic orbit becomes unstable when the absolute value of the control parameter $C$ is high enough. Eventually, for each case, the periodic orbit vanish, and the dynamics becomes chaotic. However, the road to chaos may be different according to the diffusion regime, as illustrated on Figure 4. For the cases of linear diffusion, self-diffusion on the first species, and cross-diffusion, a period doubling cascade to chaos is suspected; indeed, several period doubling are caught by the periodic orbit solver, and are visible in phase space, as a closed curve, with multiple loops, with a yellow color (indicating that orbits are revisited, it is not a transient signal). Let us note that the period doubling cascade observed in the linear case differs from the quasiperiodicity route originally observed by BVAM. This could be due to the different boundary conditions under consideration: homogeneous Dirichlet in the original case, as opposed to periodic condition in this article. For the case of self-diffusion on the second species, a fold bifurcation is observed. Observing the evolution of the density, it seems that the periodic orbit, consisting of three oscillating stripes, evolves towards a five stripes configuration, and, doing so, becomes chaotic.
<figure><img src="image_8.png" /><figcaption>Figure 4: Road to chaos, from the instability of the periodic orbit (Figure 3). $C = -1.531, -1.49, -1.63, -1.52$. A period doubling cascade appears in three cases: linear diffusion, self-diffusion on $u_1$, and cross-diffusion, clearly visible in phase space, where the closed orbits double.</figcaption></figure>

### 4.3. Strange attractors

After a transient period $\tau$, the limit set $A$ can take various forms, depending on the value of parameter $C$. Before the Hopf bifurcation, $A$ is reduced to a point. Passed the bifurcation, a periodic orbit appears, then $A$ is a circle. Passed the period doubling bifurcation, $A$ is a connected sum of two embedded circles. In the chaotic regime, $A$ fills a volume of the phase space, forming a strange attractor.

Results are displayed on Figure 5. In each of the four regimes of diffusion, we can observe similarities. The mix of blue and yellow colors, on the local phase space, is an indication of orbit revisit. Therefore orbits are bounded, and visited repeatedly through time. The power spectrum is continuous (as opposed to clear separated lines), therefore the signal is not periodic. Clear differences appear on the qualitative description of the attractor in energy space.

## 5. Conclusion

In this article, we have introduced an original model (1) which contains a rich dynamics, exhibiting original Turing patterns, oscillating patterns, and chaos. We have derived accurate numerical methods in order to approximate orbits, steady states, and periodic orbits of the latter. We also showed a convergence property (14). We have studied bifurcations of the model with respect to one control parameter, in four types of diffusion regimes: a linear case (that may be seen as a witness test), self-diffusion for each of the two species, and cross-diffusion between the species. In each of the four cases, we have observed a three steps
<figure><img src="image_9.png" /><figcaption>Figure 5: Chaos and strange attractors. Orbits of model (1), starting from unstable periodic orbit of (Figure 3), in each of the four diffusion regimes, with $C = -1.5329, -1.53, -1.63, -1.51$. For each orbit we display the local phase space (24), the spectrum of energy (3), and energy phase space (25). We observe the development of a strange attractor in every case.</figcaption></figure>
development. Initially, a stable Turing pattern appears. Then, as the control
parameter decreases, stability is lost through a Hopf bifurcation (Figure 1). A
periodic orbit appears, leading to an oscillating Turing pattern. Eventually, as
the control parameter decreases further, the periodic orbit loses its stability via
a period doubling bifurcation (Figure 2). Decreasing further the control param-
eter, the double-period orbit itself loses stability (Figure 3), but, interestingly,
scenarios differ according to different regimes. A period doubling cascade is ob-
served (Figure 4) in the linear case, as well as in the first self-difusion case and
the cross-diffusion case. To the best of author’s knowledge, it is the first time
that a clear raod to chaos is described, via bifurcation analysis, in the context
of self- and cross-diffusion. In the second self-diffusion case, a fold bifurcation
is observed. Eventually, a chaotic regime is reached in all cases, leading to the
formation of chaotic attractors, differing according to the nonlinear regime that
produced them. During this study, we have used an original energetic approach,
following [38]. The latter has revealed to be more global than the local study,
especially when studying the strange atttractors (Figure 5), allowing notably a
much clearer separation in terms of strange attractor topology.

As further works, several paths may be explored. In the regime of self-diffusion
of the second species, the road to chaos has not been clearly identified, and
therefore several options are possible, from the fold bifurcation point observed
in Figure 3. Another interesting question, that arise naturally when comparing
the strange attractors in Figure 5, would be to analyze their topology, using a
branched manifold approach, in the spirit of [36], and find their classification in
the taxonomy proposed in [37].

References

[1] Pierre-Simon de Laplace, "Essai philosophique sur les probabilités", (1813).

[2] James Clerk Maxwell, "Matter and Motion", (1876).

[3] Gaston Floquet, "Sur les équations différentielles à coefficients périodiques", Annales Ecole Normale Supérieure, **12**, 47-88, (1883).
[4] Henri Poincaré, "Les méthodes nouvelles de la Mécanique Céleste", Gauthier-Villars, Paris, (1892).

[5] Duffing, G. "Erzwungene schwingungen bei veränderlicher eigenfrequenz", F. vieweg u. Sohn, Braunschweig, (1918).

[6] B. Van der Pol et J. Van der Mark, "Frequency demultiplication", Nature, **120**, 363-364, (1927).

[7] Eberhard Hopf, "Bifurcation of a periodic solution from a stationary solution of system of differential equations", (1942).

[8] Alan Turing, "The chemical basis of morphogenesis", Philosophical Transactions of the Royal Society of London B, **237**, 37-72, (1952).

[9] Chirikov, "Resonance processes in magnetic traps", Atomic Energy, **6**, (1959).

[10] Edward Lorenz, "Deterministic nonperiodic flow", Journal of the Atmospheric Society, **20**, 130-141, (1963).

[11] Robert May, "Simple mathematical models with a very complicated dynamics", Nature, **261**, 459-467, (1976).

[12] Rössler, Otto, "An equation for chemical chaos", Physical Letters, **57 A**, 397-398, (1976).

[13] Mackey, M. C. and Glass, L., "Oscillation and chaos in physiological control systems", Science, 197(4300):287-289, (1977).

[14] Kuramoto, Y., "Diffusion induced chaos in reaction systems", Progress of Theoretical Physics Supplement, **64**, 146-367, (1978).

[15] Mitchell Feigenbaum, "Quantitative universality for a class of non-linear transformations", Journal of Statistical Physics, **19**, 25-52, (1978).

[16] Coullet, P., Tresser, C., "Itération d'endomorphismes et groupe de renormalisation", J. Phys. Colloque, C, 25-28, (1978).
[17] N. Shigesada, K. Kawasaki and E. Teramoto, "Spatial segregation of interacting species", Journal of Theoretical Biology, **79**, 83-99, (1979).

[18] Linsay, Paul S. "Period doubling and chaotic behavior in a driven anharmonic oscillator", Physical Review Letters, 47(19), 1349-1352, (1981).

[19] Libchaber, A., Laroche, C. Fauve, S. "Period doubling cascade in mercury, a quantitative measurement", Journal de Physique Lettres, 43(7), 211-216, (1982).

[20] Ermentrout, B., "Period doubling and possible chaos in neural models", SIAM Journal on Applied Mathematics, 44(1), 80-95, (1984).

[21] Papageorgiou, D.T., Smyrlis, Y.S., "The route to chaos for the Kuramoto-Sivashinsky equation", Theoretical and Computational Fluid Dynamics, 3(1), 15:42, (1991).

[22] James Yorke, "Chaos in a double pendulum", American Journal of Physics, (1991).

[23] Popov, Fioretti, "Structure of chaos in laser with saturable", Physical Review Letter, (1992).

[24] Chua, L. O., "The genesis of Chua circuit", Archiv fur elektronik und ubertragung-stechnik, (46):250-257, (1992).

[25] C. Gasquet and P. Witomski, "Fourier Analysis and Applications: Filtering, Numerical Computation, Wavelets" Springer, (1999).

[26] Barrio, R., Varea, C., Aragon, J., Maini, P., "A two dimensional numerical study of spatial pattern formation in interacting Turing systems", Bulletin of Mathematical Biology, 61(3), 483-505 (1999)

[27] I. Epstein and Vanag, "Oscillatory Turing patterns in reaction-diffusion systems with two coupled layers", Physical Review Letter, **90**, (2003).
[28] Manneville, P. [2004] "Instabilities, chaos and turbulence", *Imperial College Press*, London.

[29] Tuckerman, L. S., Huepe, C., & Brachet, M. E. [2004] "Numerical methods for bifurcation problems. In Instabilities and Nonequilibrium Structures IX", *Springer*, Netherlands.

[30] Y. Kuznetsov, "Elements of applied bifurcation theory", Springer New-York, (2004).

[31] Jean-Pierre Demailly, "Numerical analysis and differential equations", EDP Sciences, (2004).

[32] D.A. Knoll, D.E. Keyes, "Jacobian free Newton-Krylov methods: a survey of approaches and applications", Journal of Computational Physics, **193**(2), (2004).

[33] J. L. Aragon, R. A. Barrio, T. E. Woolley, R. E. Baker, and P. K. Maini, "Nonlinear effects on Turing patterns: Time oscillations and chaos", Phys. Rev. E, **86**, (2012).

[34] Pierre Collet, "A short historical account of period doublings in the pre-renormalization era", Comptes Rendus Mécanique, **347**, 287-293, (2019).

[35] Moussa, A., Perthame, B. and Salort, D., "Backward Parabolicity, Cross-Diffusion and Turing Instability", Journal of Nonlinear Science, **29**, 139-162, (2019).

[36] Gisela D. C., Guillermo, A. and Sciamarella, D., "Topology of dynamical reconstructions from Lagrangian data", Physica D: Nonlinear Phenomena, **405**, 1323-1371, (2020).

[37] Letellier, C., Stankevich, N. and Rössler, O., "Dynamical Taxonomy: Some Taxonomic Ranks to Systematically Classify Every Chaotic Attractor", International Journal of Bifurcation and Chaos, **32**, 2230004, (2022).
[38] B. Aymard, "On pattern formation in reaction-diffusion systems containing self-cross diffusion terms", Communications in Nonlinear Science and Numerical Simulation, **105**, 106090, (2022).

[39] B. Aymard, "Bifurcation analysis and steady-state patterns in reaction-diffusion systems augmented with self- and cross-diffusion", International Journal of Bifurcations and Chaos, **33**, 2350102, (2023).