# From Turing patterns to chimera states in the 2D Brusselator model

A. Provata<EMAIL_ADDRESS>Institute of Nanoscience and
Nanotechnology, National Center for Scientific Research “Demokritos”, 15341
Athens, Greece

###### Abstract

The Brusselator has been used as a prototype model for autocatalytic
reactions, and in particular for the Belouzov-Zhabotinsky reaction. When
coupled at the diffusive limit, the Brusselator undergoes a Turing bifurcation
resulting in the formation of classical Turing patterns, such as spots,
stripes and spirals in 2 spatial dimensions. In the present study we use
generic nonlocally coupled Brusselators and show that in the limit of the
coupling range $R\to 1$ (diffusive limit), the classical Turing patterns are
recovered, while for intermediate coupling ranges and appropriate parameter
values chimera states are produced. This study demonstrates how the parameters
of a typical nonlinear oscillator can be tuned so that the coupled system
passes from spatially stable Turing structures to dynamical spatiotemporal
chimera states.

††preprint: AIP/123-QED

> Oscillating reactions, also known as chemical clocks, have been studied
> intensively in the past, with notable examples the Bray-Liebhafsky reaction
> (1921), the Belousov-Zhabotinsky reaction (1958), the Briggs-Rauscher
> reaction (1972) and the Chlorite-Iodide-Malonic Acid (CIMA) reaction (1990).
> They are characterized by periodical color changing or periodic gas
> emissions and this periodic activity was earlier considered as a
> contradiction to common sense. The explanation of the mechanisms behind the
> unexpected concentration oscillations was proposed in the 1970’s, and relies
> on the nonlinear nature of the interactions between the species, which give
> rise to a limit cycle following a Hopf bifurcation scenario Epstein and
> Pojman (1998). The same nonlinear mechanisms are responsible for the
> emergence of spectacular Turing patterns when spatial diffusion is taken
> into account, as demonstrated by Alan Turing in his 1952 seminal paper
> Turing (1952). In the present study we quest the existence of a different
> unexpected spatiotemporal effect, the chimera states, which are mainly
> observed in nonlinear reactive dynamics at intermediate coupling ranges,
> whereas the Turing patterns emerge at the limit of local diffusive coupling.
> To that purpose, we study the generic Brusselator model Nicolis and
> Prigogine (1977); Lefever and Nicolis (1971), which was proposed as a
> reduction of the complex chemical mechanism representing the Belouzov-
> Zhabotinsky reaction. For this model we show, in a unifying description,
> that by appropriate choice of the parameter values and coupling ranges we
> can pass from the spatial Turing patterns to the oscillating spatiotemporal
> chimera states. In the simulations coupled Brusselators linked in 2D torus
> geometry are employed, where both Turing patterns and chimera states
> demonstrate spectacular formations.

##  I Introduction

The motivation of the present study is to describe, in a unifying way, the
emergence of spatial structures and spatiotemporal patterns in different
regions of the parameter space in the same dynamical system. For this reason
we use a typical nonlinear dynamical system with minimal number of variables
(two), the Brusselator, constructed specifically to model the complex
spatiotemporal organization of the Belouzov-Zhabotinsky reaction (see Nicolis
(1995) and references therein). For this model, we numerically demonstrate
that it is possible to drive the system from spatial Turing pattern to
spatiotemporal oscillating chimera states by controlling the range of the
diffusion from classical (local) diffusion to intermediate range (nonlocal)
diffusion with appropriate choice of the reaction rates (parameters).

Turing patterns were first predicted in 1952 by Alan Turing in a set of
nonlinear PDEs describing a hypothetical chemical reaction, which involved
autocatalytic steps and species diffusion Turing (1952). Experimentally, these
states were observed in several reactive systems beyond the Belouzov-
Zhabotinsky, such as the Bray–Liebhafsky, the Briggs–Rauscher and the
Chlorite-Iodide-Malonic Acid (CIMA) reactions Epstein and Pojman (1998);
Epstein and Showalter (1996); Mikhailov and Showalter (2006); Winfree (1978);
Kepper  _et al._ (1991); Castets  _et al._ (1990). Regarding the numerical
evidence of Turing patterns, until the late 1980’s it was generally admitted
in the literature that “Albeit their prediction by Turing, chemical periodic
stationary patterns under far from equilibrium conditions, … have not yet been
observed in a clear-cut way” Borckmans  _et al._ (1987). Nevertheless, a few
years later, in the early 1990s, numerical simulations corroborated the
experimental findings and the theoretical predictions. Turing patterns were
then reported using a number of models, including the Brusselator Verdasca
_et al._ (1992); Borckmans  _et al._ (1992); Pena and Perez-Garcia (2001), the
two- and three-variable Oregonator models Field  _et al._ (1972); Field and
Noyes (1974); Becker and Field (1985); Berenstein  _et al._ (2014); Hildebrand
_et al._ (2001); Epstein and Showalter (1996) and the two-variable
Zhabotinsky-Buchholtz-Kiyatkin-Epstein (ZBKE) model Zhabotinsky  _et al._
(1993); Nkomo  _et al._ (2013).

On the other end, chimera states have recently been the subject of intensive
studies since their first discovery by Kuramoto and Battogtokh in 2002
Kuramoto and Battogtokh (2002); Kuramoto (2002). The term “chimera”, proposed
by Abrams and Strogatz Abrams and Strogatz (2004), dominated two years after
their first introduction and refers precisely to a hybrid state. Indeed, in
networks composed of identical and identically linked oscillators, as Kuramoto
and Battogtokh demonstrated, a hybrid state is often established which
consists of stable, coexisting, coherent and incoherent domains. Originally,
the chimera states were observed under nonlocal network connectivity, but
later works using local connectivity and non-identical oscillators have also
reported emerging chimera states. In the first studies, the phase oscillator
was used as the dynamics on the network nodes Bolotov  _et al._ (2016);
Goldschmidt  _et al._ (2019); Bolotov  _et al._ (2020), while most recent
studies have demonstrated the presence of chimera states in a number of
nonlinear oscillators, such as the Hodgkin-Huxley (HH) Sakaguchi (2006), the
FitzHugh-Nagumo (FHN) Omelchenko  _et al._ (2013, 2015); Schmidt  _et al._
(2017); Zakharova  _et al._ (2017); Shepelev and Vadivasova (2019), the
Hindmarsh-Rose Hizanidis  _et al._ (2014), the Van der Pol Ulonska  _et al._
(2016), the Stuart-Landau Yeldesbay  _et al._ (2014); Gjurchinovski  _et al._
(2017); Tumash  _et al._ (2018), the Leaky Integrate-and-Fire (LIF) Olmi  _et
al._ (2015); Tsigkri-DeSmedt  _et al._ (2017, 2018) and other oscillator
networks. In the above studies the presence of chimera (multichimera) states
were reported with one (many) coherent and incoherent domains. The patterns
and multiplicity of the chimeras depend on the parameters governing both the
nodal dynamics and the network coupling.

The role that the network architecture plays in the form of the chimera states
has also been the subject of recent intense investigations. Most of the
studies cited above where performed in a 1D connectivity with periodic
boundary condition (K-ring connectivity). Considerably fewer studies concern
connectivities in higher spatial dimensions. Studies in 2D-square (torus)
connectivity have revealed coherent (and incoherent) domains of circular disc
shapes, stripes, spots, scroll waves etc Schmidt  _et al._ (2017); Maistrenko
_et al._ (2017); Laing (2017); Omel’Chenko and Knobloch (2019). Similar
studies in 3D-cubic (hypertorus) connectivity produce coherent and incoherent
domains of higher complexity, such as spheres, cylinders, layers and
combinations of these Omel’Chenko and Knobloch (2019); Maistrenko  _et al._
(2015); Kasimatis  _et al._ (2018); Maistrenko  _et al._ (2017, 2020);
Koulierakis  _et al._ (2020).

Almost one decade after the first theoretical prediction of chimera states,
experimental evidence of these states were reported in a number of mechanical,
physical and chemical systems consisting of interacting oscillatory units.
Some examples reported in the literature include chimeras in experiments in
optical systems M.Hagerstrom  _et al._ (2012), in electronic circuits Gambuzza
_et al._ (2014), in mechanics (systems of coupled mechanical oscillators)
Martens  _et al._ (2013); Kapitaniak and Kurths (2014); Dudkowski  _et al._
(2016); Manoj  _et al._ (2021), in biomedicine Cherry and Fenton (2008);
Mormann  _et al._ (2000, 2003); Santos  _et al._ (2017), and in reaction
diffusion systems Tinsley  _et al._ (2012); Wickramasinghe and Kiss (2013);
Schmidt  _et al._ (2014). Beyond experimental findings in the laboratory,
chimera states have been associated with the uni-hemispheric sleep in mammals
and birds Rattenborg  _et al._ (2000); Rattenborg (2006); Ramlow  _et al._
(2019), with the onset of epileptic seizures Mormann  _et al._ (2000, 2003)
and other biomedical conditions Andrzejak  _et al._ (2016); Ruzzene  _et al._
(2019); Bansal  _et al._ (2019); Koulierakis  _et al._ (2020). For reviews on
the recent theoretical and experimental developments in the field of chimera
states, see references Panaggio and Abrams (2015); Schöll (2016); Yao and
Zheng (2016); Majhi  _et al._ (2019); Omel’Chenko (2018); Zakharova (2020);
Omel’Chenko (2022).

Regarding the emergence of chimera states in reaction-diffusion systems, the
first such reports date as early as 2012. In particular, photochemically
coupled Belousov-Zhabotinsky micro-oscillators have been shown to produce
complex synchronization phenomena, such as spiral wave chimeras, spiral wave
core splitting and phase cluster states Tinsley  _et al._ (2012); Nkomo  _et
al._ (2013); Taylor  _et al._ (2015); Nkomo  _et al._ (2016); Totz  _et al._
(2018, 2020). Numerical simulations using the dimensionless Zhabotinsky-
Buchholtz-Kiyatkin-Epstein (ZBKE) model with a time delay in the coupling
between the oscillators corroborate the experimental findings Nkomo  _et al._
(2016); Totz  _et al._ (2018, 2020). These studies indicate that chimera
states are not only theoretical predictions but are realized in chemical
experimental setups. Additional evidence of chimera states in reaction-
diffusion systems has been reported by the present author and collaborators in
a different numerical study drawing from population dynamics and reactive
systems in cyclic reactions Hizanidis  _et al._ (2015). In this study, the
Lattice Limit Cycle (LLC) model was considered with 4-th order nonlinearities.
Nonlinearly coupled LLCs were shown to produce chimera states in the spiking
regime, when the bifurcation parameter increases away from the Hopf
bifurcation. More recent studies investigating the presence of Turing patterns
and chimera states in networks of nonlinear oscillators have been reported in
the literature for the phase oscillator Franović  _et al._ (2021) and for
coupled SQUID models Hizanidis  _et al._ (2020).

As earlier stated, in the present study we consider a well known reaction
diffusion scheme, the Brusselator, proposed as a minimal prototype scheme for
demonstrating the effects of nonlinear interactions in chemistry. Its minimal
dynamics involves two interacting variables with 3rd order nonlinearities and
an autocatalytic step. Numerous studies in the past have shown that this model
gives rise to interesting Turing patterns in 1D, 2D and 3D spatial geometries,
when normal diffusion in the form of nearest neighbor interactions is included
in the dynamics and the parameter values are chosen in the Turing bifurcation
regime Borckmans  _et al._ (1992); Verdasca  _et al._ (1992); Pena and Perez-
Garcia (2001); Kokolnikov  _et al._ (2006); Budroni and Wit (2016); Bullara
_et al._ (2013); Tildi  _et al._ (2016); Shoji  _et al._ (2007). In the
sequel, we extend gradually the range of interactions, going from nearest
neighbors (normal diffusion) to nonlocal, intermediate ranges where the
chimera states prevail. In parallel, we determine the reaction rates and
parameter ranges which trigger the appearance of complex synchronization
phenomena and chimera states in two spatial dimensions.

The organization of the study is as follows: In the next section, II, we
introduce the uncoupled (II.1) and the coupled (II.2) Brusselator model and
discuss its mean field dynamics and the conditions and parameters for the Hopf
and the Turing bifurcations. At the end of the same section we introduce some
quantitative measures appropriate for the description of Turing patterns and
chimera states. In section III, we show that by taking the diffusing limit in
the coupling range we obtain Turing structures of different morphologies. In
section IV, we increase the coupling range and for appropriate parameter
values (after Hopf bifurcation) we demonstrate the emergence of chimera
states. By scanning the coupling range and coupling angle parameters we
determine the region in the parameter space where the chimera states dominate.
In section V, the nature of Turing and chimera patterns are discussed using
the Fourier transforms of the simulated evolution of the systems’ variables.
In the concluding section, we recapitulate our main results and propose open
problems.

##  II The Model

In this section we first present the ODEs of the uncoupled Brusselator model
and discuss, in brief, the Hopf bifurcation scenario giving rise to a limit
cycle. Next, we present the nonlocally coupled Brusselator network and discuss
the emergence of a Turing bifurcation at the limit of local coupling. In the
last part of this section, we introduce some quantitative indices for
quantifying the spatial and temporal properties of the system.

###  II.1 The uncoupled Brusselator Model

As discussed in the Introduction, the Brusselator model was introduced in late
1960’s by R. Lefever, G. Nicolis and I. Prigogine as a minimal, two-variable,
autocatalytic model producing all spatiotemporal complexity of the Belouzov-
Zhabotinsky oscillatory reaction Nicolis and Prigogine (1977); Lefever and
Nicolis (1971). It consists of two ODEs of 3rd order, with two variables,
$\displaystyle X$ and $\displaystyle Y$, and two control parameters,
$\displaystyle A$ and $\displaystyle B$, as follows:

$\displaystyle\displaystyle\frac{dX}{dt}$
$\displaystyle\displaystyle=A-(B+1)X+X^{2}Y$ (1a)
$\displaystyle\displaystyle\frac{dY}{dt}$
$\displaystyle\displaystyle=BX-X^{2}Y.$ (1b)

In Eqs. (1), the variables $\displaystyle X$ and $\displaystyle Y$ take
positive values since they were first constructed to represent chemical
species concentrations, while the parameters $\displaystyle A$ and
$\displaystyle B$ represent the initial product concentrations and remain
constant (considered as sources).

The Brusselator model, Eqs. (1), has a single fixed point, $\displaystyle
X_{s}$,

$\displaystyle\displaystyle(X_{s},Y_{s})=(A,B/A).$ (2)

Stability analysis around $\displaystyle(X_{s},Y_{s})$ reveals that this fixed
point changes its asymptotic stability when the parameter $\displaystyle B$
attains its critical value $\displaystyle B_{c}$:

$\displaystyle\displaystyle B_{c}=A^{2}+1.$ (3)

Namely, for $\displaystyle B<B_{c}=A^{2}+1$, $\displaystyle(X_{s},Y_{s})$ is a
stable node or stable focus and the phase space trajectories are attracted to
the fixed point. At $\displaystyle B=B_{c}$ the dynamics around the fixed
point changes from asymptotic stability to instability and beyond this point a
limit cycle is created Nicolis (1995). Because limit cycles are stable
oscillators under reasonable perturbations, they are prerequisites for the
emergence of the complex synchronized states, chimera states, when many limit
cycles are coupled in networks. For this reason, in section IV the parameter
values $\displaystyle A=1.0$ and $\displaystyle B=2.5-3.5$ will be used (such
that $\displaystyle B>A^{2}+1$), when the network composed by Brusselator
units will be tuned to produce chimera states.

On the other hand, in the Brusselator model the Turing structures are produced
away from the regime of limit cycles and under the Turing bifurcation
conditions, where the coupling terms need to be taken into account. For this
reason, the discussion on the precise Brusselator parameter values which lead
to the emergence of Turing structures is postponed till section III, after the
coupling terms have been included.

###  II.2 The coupled Brusselator Model

In this section, spatial coupling is introduced between units in 2D, square
lattice geometry. Earlier studies have discovered spectacular Turing patterns
in 2D, such as dots and stripes, and the aim of the present work is to clarify
the different parameter domains in the 2D Brusselator network where the Turing
patterns and the chimera states reside. On the square lattice network
containing $\displaystyle N\times N$ units, the Brusselator located at
position coordinates $\displaystyle(i,j)$ is described by the time dependent
variables $\displaystyle X_{ij}$ and $\displaystyle Y_{ij}$. The communication
between the Brusselators is achieved by the coupling matrix $\displaystyle C$
whose elements, $\displaystyle c_{ij:kl}$, represent the strength of the link
between Brusselators at positions $\displaystyle(i,j)$ and
$\displaystyle(k,l)$ on the network, with $\displaystyle i,j,k,l\leq N$.
Recent investigations of synchronization phenomena on complex networks have
shown that coupling matrices with strong cross-coupling terms
$\displaystyle(w_{XY},\>w_{YX})$ and both excitatory (positive) and inhibitory
(negative) terms are necessary (but not sufficient) for partial
synchronization in the form of chimera states Omelchenko  _et al._ (2013);
Omel’Chenko (2018). Under these conditions, the nonlocally coupled Brusselator
dynamics on the 2D square lattice is governed by the following equations:

$\displaystyle\displaystyle\frac{dX_{ij}}{dt}$
$\displaystyle\displaystyle=A-(B+1)X_{ij}+X_{ij}^{2}Y_{ij}+\sum_{(k,l)}c_{kl:ij}\left[w_{XX}(X_{kl}-X_{ij})+w_{XY}(Y_{kl}-Y_{ij})\right]$
(4a) $\displaystyle\displaystyle\frac{dY_{ij}}{dt}$
$\displaystyle\displaystyle=BX_{ij}-X_{ij}^{2}Y_{ij}+\sum_{(k,l)}c_{kl:ij}\left[w_{YX}(X_{kl}-X_{ij})+w_{YY}(Y_{kl}-Y_{ij})\right]$
(4b)

In Eqs. 4, we implicitly assume that all Brusselators are identical since they
have the same parameters $\displaystyle A$ and $\displaystyle B$.

The coupling matrix $\displaystyle C$ is limited by the range $\displaystyle
R$ and strength $\displaystyle\sigma$ of the nonlocal interactions. The
elements $\displaystyle c_{ij:kl}$ representing the interactions between
Brusselators at positions $\displaystyle(i,j)$ and $\displaystyle(k,l)$ have
the form:

$\displaystyle\displaystyle
c_{ij:kl}=\begin{cases}\dfrac{\sigma}{\displaystyle{(2R+1)^{2}-1}},&\text{when}\>\>i-R\leq
k\leq i+R,\>\>\>j-R\leq l\leq j+R\\\
\>\>\>\>\>\>\>\>\>\>\>\>\>\>0,&\text{otherwise}.\end{cases}$ (5)

Namely, the Brusselator at position $\displaystyle(i,j)$ is coupled with all
other elements located only within a rectangle of size
$\displaystyle(2R+1)\>{\rm X}\>(2R+1)$, whose center is at position
$\displaystyle(i,j)$. Each nonzero link carries constant weight (coupling
strength) $\displaystyle\sigma$, while the overall contribution is normalized
by the total number of nonzero links, which is $\displaystyle(2R+1)^{2}-1$. In
Eqs. (4) and (5), periodic boundary conditions on the torus are considered
and, therefore, all indices $\displaystyle i,j,k,l,$ are taken
$\displaystyle\mod N$.

Apart for the coupling matrix $\displaystyle C$ between the nodes of the
network, a rotational matrix, $\displaystyle W$, is introduced to account for
the cross-coupling terms. $\displaystyle W$ contains both diagonal coupling,
in the form $\displaystyle(w_{XX},w_{YY})$ and cross coupling
$\displaystyle(w_{XY},w_{YX})$ and is parameterized by a single parameter
$\displaystyle\varphi$, which is close to $\displaystyle\pi/2$ if cross-
coupling is dominant Omelchenko  _et al._ (2013). This coupling phase
$\displaystyle\varphi$ is similar to the phase-lag parameter
$\displaystyle\alpha$ of the paradigmatic Kuramoto phase oscillator model,
which is widely used to generically describe coupled oscillator networks. The
coupling phase is necessary for the development of nontrivial partial
synchronization patterns, as has been shown for the Kuramoto model Omel’Chenko
_et al._ (2010) and for the FHN model Omelchenko  _et al._ (2013). This
rotational matrix $\displaystyle W$ is also used in Refs. Omelchenko  _et al._
(2015); Schmidt  _et al._ (2017); Ramlow  _et al._ (2019); Mikhaylenko  _et
al._ (2019) and has the form:

$W=\begin{pmatrix}w_{XX}&w_{XY}\\\
w_{YX}&w_{YY}\end{pmatrix}=\begin{pmatrix}\cos{\varphi}&\sin{\varphi}\\\
-\sin{\varphi}&\cos{\varphi}\end{pmatrix}.$ (6)

In the present study different values of the coupling range $\displaystyle R$
and the coupling phase $\displaystyle 0\leq\varphi\leq 2\pi$ are examined to
determine the domains in the $\displaystyle(\varphi,R)$ parameter intervals
where the Turing patterns or the chimera states dominate.

###  II.3 Quantitative measures

The synchrony/asynchrony in the coupled Brusselator system can be verified
using the mean phase velocity profile, as in most studies where many
oscillators are interlinked. The mean phase velocity
$\displaystyle\omega_{ij}$ accounts for the average number of cycles that
oscillator at position $\displaystyle(i,j)$ has performed in the time unit
Omelchenko  _et al._ (2015); Tsigkri-DeSmedt  _et al._ (2018). Let us call
$\displaystyle Q_{ij}$ the number of complete cycles that
oscillator$\displaystyle(i,j)$ has performed during a certain time interval
$\displaystyle\Delta T$. Then $\displaystyle\omega_{ij}$ is defined as:

$\displaystyle\displaystyle\omega_{ij}=2\pi\frac{Q_{ij}}{\Delta T}=2\pi
f_{ij}\>\>,$ (7)

where $\displaystyle f_{ij}$ is the frequency of the oscillator at position
$\displaystyle(i,j)$. Because the definitions of frequency and mean phase
velocity differ by just a factor of $\displaystyle 2\pi$, in the following the
two expressions (“frequency” and “mean phase velocity”) are used
interchangeably.

Chimera states are characterized by a distribution of frequencies, while the
coherent states are characterized by a single, common frequency. Full
incoherence may also be characterized by a single, common frequency, in some
cases. Let us denote by $\displaystyle\omega_{\rm min}=\min\\{\omega_{ij}\\}$,
$\displaystyle 1\leq i,j\leq N$, the minimum mean phase velocity recorded in
the Brusselator network. Similarly, $\displaystyle\omega_{\rm
max}=\max\\{\omega_{ij}\\}$, $\displaystyle 1\leq i,j\leq N$, is the maximum
mean phase velocity observed. The difference

$\displaystyle\displaystyle\Delta\omega=\omega_{\rm max}-\omega_{\rm min},$
(8)

is a measure of the extent of the mean phase velocity distribution in the
system. We may assume that if $\displaystyle\Delta\omega>2*\pi/\Delta T$, then
the difference of the mean phase velocity between nodes can be considered
significant and may be an indication of the establishment of a chimera state
in the network. If $\displaystyle\Delta\omega<2*\pi/\Delta T$, this might be
considered as a one-cycle over- or under-estimation and can be regarded as the
error threshold.

The Kuramoto order parameter is another measure often used to quantify
coherence/incoherence in the system Kuramoto and Battogtokh (2002); Bick  _et
al._ (2016). Let us denote $\displaystyle\theta_{ij}$ the instantaneous phase
of Brusselator at position $\displaystyle(i,j)$, which is defined as
$\displaystyle\theta_{ij}=\arctan{\left[Y_{ij}/X_{ij}\right]}$. The Kuramoto
order parameter $\displaystyle Z$ is defined globally, over the entire system
as:

$Z=\frac{1}{N^{2}}\displaystyle{\left|\sum_{i,j}^{N}\mathrm{e}^{+i\theta_{ij}}\right|}.$
(9)

In Eq. (9), if the system is synchronous and all phases are identical, then
$\displaystyle Z=1$. If all phases are random as in the asynchronous state,
the various phases in the sum of Eq. (9) cancel out and $\displaystyle Z=0$.
For intermediate cases, where synchronous and asynchronous regions are
present, $\displaystyle 0<Z<1$, and this is the domain where chimera states
can be observed. On the other hand, even in cases that the oscillators in the
network have identical mean phase velocities but they operate under constant
phase lags, they can be mistakenly considered as chimera states, because
$\displaystyle Z<1$. Therefore, both measures $\displaystyle\Delta\omega>0$
and $\displaystyle 0<Z<1$ are necessary to quantify synchronization in
networks of coupled elements.

Other measures such as the spatial and/or temporal Fourier spectra can be used
to further verify the spatial and/or temporal order of the network nodes,
respectively.

##  III Turing Structures in 2D

Classical Turing structures are obtained when a nonlinear reactive system (the
Brusselator in the present case) is subjected to interactions in the form of
diffusion. The diffusion terms which are added to the original Brusselator
model, Eqs. (1), involve the Laplacian operator and have the form
$\displaystyle D_{X}\nabla^{2}X$ for Eq. (1a) and $\displaystyle
D_{Y}\nabla^{2}Y$ for Eq. (1b). The parameters $\displaystyle D_{X}$ and
$\displaystyle D_{Y}$ represent the diffusion coefficients and Eqs. (1) then
become:

$\displaystyle\displaystyle\frac{dX}{dt}$
$\displaystyle\displaystyle=A-(B+1)X+X^{2}Y+D_{X}\nabla^{2}X$ (10a)
$\displaystyle\displaystyle\frac{dY}{dt}$
$\displaystyle\displaystyle=BX-X^{2}Y+D_{Y}\nabla^{2}Y.$ (10b)

The study of these equations demonstrates that the fixed point
$\displaystyle(X_{s},Y_{s})$ undergoes a Turing bifurcation at a critical
value $\displaystyle B^{T}_{c}$ of the parameter B, such thatNicolis (1995);
Verdasca  _et al._ (1992); Pena and Perez-Garcia (2001):

$\displaystyle\displaystyle B^{T}_{c}$ $\displaystyle\displaystyle=$
$\displaystyle\displaystyle(1+\eta A)^{2},$ $\displaystyle\displaystyle\eta$
$\displaystyle\displaystyle=$
$\displaystyle\displaystyle\left(\frac{D_{X}}{D_{Y}}\right)^{1/2}<(1+A^{2})^{1/2}-1.$
(11)

In addition, the discretization in space of the Laplacian operator retains
only interactions with the first neighbors, and Eqs. (10) reduce to (see
Appendix A for the derivation):

$\displaystyle\displaystyle\frac{dX_{ij}}{dt}$
$\displaystyle\displaystyle=A-(B+1)X_{ij}+X_{ij}^{2}Y_{ij}+\frac{D_{X}}{h^{2}}\sum_{\\{k,l\\}}\left[X_{kl}-X_{ij}\right]$
(12a) $\displaystyle\displaystyle\frac{dY_{ij}}{dt}$
$\displaystyle\displaystyle=BX_{ij}-X_{ij}^{2}Y_{ij}+\frac{D_{Y}}{h^{2}}\sum_{\\{k,l\\}}\left[Y_{kl}-2Y_{ij}\right].$
(12b)

Here the indices vary as $\displaystyle i-1\leq k\leq i+1$ and $\displaystyle
j-1\leq l\leq j+1$, (i.e., the coupling range $\displaystyle R=1$) and
$\displaystyle h$ denotes the lattice constant. Equations (12) are of the same
form as Eqs. (4), with the following parameter identifications

  * a)

$\displaystyle c_{ij;kl}\>w_{XX}=D_{X}/h^{2}$,

  * b)

$\displaystyle c_{ij;kl}\>w_{YY}=D_{Y}/h^{2}$,

  * c)

$\displaystyle w_{XY}=w_{YX}=0$ and

  * d)

$\displaystyle R=1$.

With these parameter identifications and under restrictions Eq. (11), the
generic coupled Brusselator system, Eqs. (4), may produce Turing patterns as
shown in Fig. 1.

For the simulations in Fig. 1 all systems started with the same random initial
conditions. As working parameter set, generic values were used (obeying
conditions Eq. (11)), namely, $\displaystyle A=4.5$, $\displaystyle
D_{X}=7.0$, $\displaystyle D_{Y}=56.0$, $\displaystyle w_{XY}=w_{YX}=0$. The
constant $\displaystyle h=1$ as it is representing the lattice constant
(absolute distance between successive nodes), as described in the Appendix A.
On the left panels of Fig. 1 the parameter $\displaystyle B$ takes the value
$\displaystyle B=6.75$, while on the right ones $\displaystyle B=8.0$. All
systems were integrated for 1500 Time Units (TUs), i.e., for $\displaystyle
1500\times 10^{3}$ integration steps. The Turing patterns stabilized between
300 and 400 TUs and remained stable thereafter.

Different values of coupling range are used, $\displaystyle R=1$ at the top
row, $\displaystyle R=2$ in the middle row and $\displaystyle R=3$ in the last
row, to demonstrate how the Turing patterns dissolve as the coupling range
increases from local (nearest neighbor) coupling to nonlocal coupling. Note
that the local coupling $\displaystyle R=1$ corresponds directly to the
diffusion operator, while for $\displaystyle R>1$ we step away from the
dominance of the diffusion coefficient and the Turing bifurcation analysis is
no longer valid.

In Fig. 1 panel (a1), stable Turing spot patterns are recorded for
$\displaystyle B=6.75$ and $\displaystyle R=1$ and all other parameters as in
the working parameter set, described in the previous paragraph and also
provided in the legend of Fig. 1. Note that the size of spots increases with
$\displaystyle D_{X}$ and $\displaystyle D_{Y}$ (not shown). In panel (a2)
stable stripe patterns are recorded, for $\displaystyle B=8.0$ and
$\displaystyle R=1$. In both panels (a1) and (a2) the parameters obey
conditions Eqs. (11) while the coupling range $\displaystyle R=1$, attains the
limit of diffusion on lattice Verdasca  _et al._ (1992); Borckmans  _et al._
(1992); Pena and Perez-Garcia (2001).

The criterion on the coupling range is relaxed in the middle row of Fig. 1.
Here the coupling range is increased to the value $\displaystyle R=2$ and now
the size of the spots and stripes increases leading to fewer spots of larger
radius in Fig. 1(b1) and fewer stripes of larger width in Fig. 1(b2).
Increasing further the coupling range to $\displaystyle R=3$ in the last row
of Fig. 1, we note that the size of the spots increases further and they merge
into strips, see Fig. 1(c1). Note also that the range of the $\displaystyle
X-$ variable decreases inversely with $\displaystyle R$ and the system tends
to the pure homogeneous state as $\displaystyle R$ increases. Similarly, in
Fig. 1(c2) the width of the stripes increases with $\displaystyle R$ and even
fewer stripes are formed, leading again to the homogeneous state as
$\displaystyle R$ increases further.

The reason for increasing of the size of the patterns (spots, stripes) with
the coupling range $\displaystyle R$ can be intuitively understood, because
for larger $\displaystyle R$ each node connects to a larger number of
neighbors and therefore the elements can adjust their variables at larger
distances/areas. A similar effect can be also seen in the next section, where
the chimera states synchronize at larger distances/areas, as the coupling
range increases.

Figure 1:  Color-coded stable snapshots ($\displaystyle X-$variable profiles),
Turing patterns, of the coupled Brusselator network. a1) Multiple spots for
$\displaystyle B=6.75$ and $\displaystyle R=1$, a2) multiple stripes, for
$\displaystyle B=8$ and $\displaystyle R=1$, b1) Multiple spots for
$\displaystyle B=6.75$ and $\displaystyle R=2$, b2) multiple stripes, for
$\displaystyle B=8$ and $\displaystyle R=2$, c1) Multiple spots for
$\displaystyle B=6.75$ and $\displaystyle R=3$ and c2) multiple stripes, for
$\displaystyle B=8$ and $\displaystyle R=3$. For all integrations Eqs. (4) and
(5) were used, with the following parameters: $\displaystyle N=64$,
$\displaystyle A=4.5$, $\displaystyle D_{X}=7.0$, $\displaystyle D_{Y}=56$,
$\displaystyle w_{XY}=w_{YX}=0$, $\displaystyle h=1$, $\displaystyle\sigma=1$
with the integration constant $\displaystyle dt=0.001$. Images are recorded at
640 TUs. Simulations start from random initial conditions with periodic
boundary conditions. Two related videos are presented in the Supplementary
Material (Multimedia view corresponding to images Fig. 1(a1)
(movie_fig1a1_aa4.5_bb6.75_du7.0_dv56.0_dt0.001.mp4) and Fig. 1(a2)
(movie_fig1a2_aa4.5_bb8.00_du7.0_dv56.0_dt0.001.mp4)).

Specific features of the patterns in Fig. 1 can be extracted via the Fourier
transform, and the spatial wave length can be considered as a function of the
coupling radius. The results are shown in Table 1.

| Spots | Stripes
---|---|---
R | Amplitude | Wavelength | Amplitude | Wavelength
1 | 1.03 | 15.84 | 2.7395(h) | 12.66(h)
2 | 0.32 | 31.79 | 2.83(h) | 21.07(h)
3 | 0.01 | 31.56 | 2.88(v) | 31.62(v)
4 | 0.58 | 63.01 | 2.97(h) | 63.10(h)
5 | $\displaystyle 7\times 10^{-8}$ | 63.01 | 2.44(h) | 63.10(h)
Table 1: The amplitude and wavelength of the spots and stripes patterns in
Fig. 1 for different values of the coupling range $\displaystyle R$. The
values related to the spot pattern are averages over the $\displaystyle x-$
and $\displaystyle y-$ spatial directions. The notation (h) or (v) in the
stripes pattern denote that the amplitude and wavelength were computed in the
horizontal or vertical direction, respectively. The parameter values are as in
Fig. 1.

In Table 1, on the 2nd and 3rd columns, the amplitudes and wavelengths of the
spot pattern are presented for different $\displaystyle R$ values, as
corresponding to Fig. 1 (left column). They were computed as averages over the
amplitudes and wavelengths in the $\displaystyle x-$ and $\displaystyle y-$
directions. It is easy to understand that as $\displaystyle R$ increases the
wavelength increases and for $\displaystyle R=4$ and 5 one wavelength covers
the entire system. In addition, for $\displaystyle R=4$ there is a non-
negligible mode amplitude ($\displaystyle=0.58$), while for $\displaystyle
R=5$ the amplitude tends to zero, indicating that the system has reached the
homogeneous steady state.

Similar conclusions are reached for the stripe pattern. Depending on the
orientation of the stripes, the amplitudes and wavelengths recorded in the 4th
and 5th columns of the Table 1 were computed either in the horizontal
($\displaystyle x-$) direction and are denoted by the letter (h) or in the
vertical ($\displaystyle y-$) direction and are denoted by the letter (v).
Here also, as $\displaystyle R$ increases from 1 to 2 and to 3, the wavelength
doubles each time. For $\displaystyle R=4$, one wavelength covers the entire
system with non-negligible amplitude ($\displaystyle=2.97$). As the coupling
range increases further to $\displaystyle R=5$, still one wavelength covers
the entire system, but the amplitude of the spatial oscillations decreases
($\displaystyle=2.44$). Therefore, the amplitude of the patterns decreases
with $\displaystyle R$, leading again to the homogeneous steady state for
large $\displaystyle R-$values.

The results of Table 1 and Fig. 1 highlight the special characteristics of the
Turing patterns. They are best observed at the diffusion limit, which
corresponds to $\displaystyle R=1$ in the discrete lattice version of the
equations, and they dissipate as $\displaystyle R$ increases leading to
homogeneous steady states. As we will see in the next section, this is not the
case for the chimera states. Not only the chimeras are observed in the
parameter regions where limit cycles are supported but they also appear away
from the diffusing limit, at intermediate coupling ranges.

We would like to stress here that we do not go into a deep study of the
pattern formation of the Brusselator model, this has been done extensively
elsewhere Borckmans  _et al._ (1992); Verdasca  _et al._ (1992); Pena and
Perez-Garcia (2001); Kokolnikov  _et al._ (2006); Budroni and Wit (2016);
Bullara  _et al._ (2013); Tildi  _et al._ (2016); Shoji  _et al._ (2007). We
only highlight the dissolution of the Turing patterns as the coupling range
increases, in contradistinction with the behavior of chimera states, which are
not supported at the diffusing limit, but only at intermediate coupling ranges
(provided that all oscillators are identical).

##  IV Chimera states in 2D

Figure 2:  Color-coded values of a) the Kuramoto order parameter
$\displaystyle Z$ and b) the frequency deviation in the system,
$\displaystyle\Delta\omega$, as a function of the coupling range
$\displaystyle R$ and the coupling angle $\displaystyle\varphi$. Other
parameters for the Brusselator model are $\displaystyle A=1.0$, $\displaystyle
B=2.9$. Integration parameters are set to $\displaystyle dt=0.010$,
$\displaystyle D_{X}=D_{Y}=1$. The network is a 2D square lattice with
periodic boundary conditions and size $\displaystyle N\times N=64\times 64$.
All simulations start from the same random initial conditions.

The same model, Eqs. (4) with connectivities as in (5) and (6), for
appropriate parameter values and intermediate coupling ranges can lead to
oscillatory steady states with coexisting variable frequencies on the network,
known as chimera states. To demonstrate this, we integrate numerically the
afore mentioned equations, using the following working parameter set:
$\displaystyle A=1.0$, $\displaystyle B=2.9$ to place the uncoupled
Brusselator in the parameter domain where limit cycles are observed, see
section II.1. The integration parameters are set to $\displaystyle dt=0.010$,
$\displaystyle D_{X}=D_{Y}=1$ and the 2D network consists of a square lattice
composed by $\displaystyle N\times N=64\times 64$ elements, with periodic
boundary conditions (same as in the previous section, III, where the Turing
patterns are produced). The coupling strength $\displaystyle\sigma$ is set to
unity. To explore the presence of chimera states in the Brusselator network,
the coupling range $\displaystyle R$ and coupling angle $\displaystyle\varphi$
were varied in the ranges, $\displaystyle 1\leq R\leq 21$ and $\displaystyle
0\leq\varphi<2\pi$ (or $\displaystyle 0\leq\varphi<360$ in degrees).

As discussed in section II.3, the presence of chimera states is characterized
by values of the Kuramoto order parameter $\displaystyle Z$ such that
$\displaystyle 0<Z<1$ and $\displaystyle\Delta\omega>0$. In Fig. 2a, we
present the values of $\displaystyle Z$ in a color-coded map, as a function of
$\displaystyle R$ and $\displaystyle\varphi$. We note that chimera states can
be sought in the regions that the $\displaystyle Z$-values take dark colors
and away from the yellow regions. Similarly, in Fig. 2b, we present the values
of $\displaystyle\Delta\omega$ in a color-coded map. In the
$\displaystyle\Delta\omega$ map, there is evidence of chimera presence in the
red and green regions which support high divergence of local frequencies in
the system.

Both images, Fig. 2a and b, give compatible information as far as the
parameter regions where chimera states are realizable. Namely, for small
values of $\displaystyle\varphi<~{}90^{\circ}$ and all values of
$\displaystyle R$ chimera states are not possible: the Kuramoto order
parameter $\displaystyle Z\sim 1$, the frequency distribution
$\displaystyle\Delta\omega\sim 0$ then coherence reins in the system.
Similarly, chimera states are not promoted for high values of
$\displaystyle\varphi>280^{\circ}$. This is in agreement with previous studies
of other models (such as the FitzHugh Nagumo model) which study and report
chimera states for $\displaystyle\varphi$ values away from 0 or $\displaystyle
2\pi$.Omelchenko  _et al._ (2013, 2015). Although there is a slight dependence
on the coupling range $\displaystyle R$, as both indices $\displaystyle Z$ and
$\displaystyle\Delta\omega$ indicate, this dependence is not very strong.
Comparing the two maps, the $\displaystyle\Delta\omega$ one (Fig. 2b) shows
higher details in the structure and higher ramifications in the frequency
distributions.

Our numerical investigations demonstrate that interesting chimera states are
produced when the values of the coherence indices escape the
$\displaystyle(Z,\Delta\omega)\sim(1,0)$ values. In particular, at the
borderlines between the yellow (black-grey) and the dark regions (red-green)
in Fig. 2a (b), one may observe the formation of stripe chimeras, several
examples of which are provided in Fig. 3. In this figure, where
$\displaystyle\varphi\approx 90^{\circ}$ and different values of
$\displaystyle R$ are used, alternation of coherent and incoherent stripes is
recorded. Judging from the color-coded maps, Figs. 2a,b, we note that these
chimera patterns are displayed for borderline parameter values, where
$\displaystyle\Delta\omega$ just exceeds 0 and $\displaystyle Z$ just falls
below 1. In Figs. 3a and c two synchronous and 2 asynchronous stripes develop,
while in panels b and d only one synchronous and one asynchronous stripe are
recorded. It is worth to note that in Fig. 3a, b, c the stripes develop
horizontally, while in d the stripes are vertical. This is due to the
evolution of the random initial conditions, which are common in all
simulations. Other pairs with notable stripe chimeras are recorded for
$\displaystyle(\varphi,R)$ values such as $\displaystyle(91^{\circ},8)$,
$\displaystyle(92^{\circ},9)$, $\displaystyle(92^{\circ},10)$,
$\displaystyle(94^{\circ},14)$, $\displaystyle(94^{\circ},15)$,
$\displaystyle(95^{\circ},15)$, $\displaystyle(96^{\circ},15)$,
$\displaystyle(95^{\circ},16)$, $\displaystyle(96^{\circ},16)$,
$\displaystyle(95^{\circ},17)$, $\displaystyle(96^{\circ},17)$ and others.

Figure 3:  Synchronous and asynchronous stripes for low
$\displaystyle\varphi$-values. Spatial $\displaystyle X-$profiles for: a)
$\displaystyle(\varphi=91^{\circ},\>\>R=9)$, b)
$\displaystyle(\varphi=93^{\circ},\>\>R=14)$, c)
$\displaystyle(\varphi=96^{\circ},\>\>R=18)$ and d)
$\displaystyle(\varphi=97^{\circ},\>\>R=19)$. Other parameters are
$\displaystyle A=1.0$, $\displaystyle B=2.9$, $\displaystyle dt=0.010$,
$\displaystyle D_{X}=D_{Y}=1$, $\displaystyle N=64$ and
$\displaystyle\sigma=1$. Toroidal periodic boundary conditions are applied.
The images are recorded at 1000 TUs. All simulations start from the same
random initial conditions. A related video is presented in the Supplementary
Material (Multimedia view corresponding to image Fig. 3c
(movie_fig03c_brusselator_sig_1.00_phi_096_R_18.mp4)).

One significant difference can be stressed between images in Fig. 1 (Turing
patterns) and those in Fig. 3 (chimera states) and the ones that follow. In
Fig. 1 the images are stationary (still) after transient and do not change
with time, while in Fig. 3 (and lateron) the images are dynamic: the colors
(state variables) oscillate in time, but the synchronous areas remain
synchronous and so do the asynchronous ones. This means that in the color-
coded spacetime plots in Fig. 3, the colors of all nodes (pixels) change
continuously between black and yellow. To the contrary, in Fig. 1 they remain
fixed after an initial transient period.

Stripes are also recorded when the high values of $\displaystyle\Delta\omega$
(low values of $\displaystyle Z$) cross to lower values of
$\displaystyle\Delta\omega$ (high values of $\displaystyle Z$); this takes
place around $\displaystyle\varphi\approx 300^{\circ}$, for different
$\displaystyle R$-values. Spatial portraits of typical stripe chimeras in
these parameter region are reported in Fig. 4. Panel a shows coexistence of a
coherent (blue) and an incoherent stripe. In panel b, the coherent stripe has
variable width. Such states are often unstable and they result in a turbulent
incoherent state. Panels c (horizontal stripes) and d (vertical stripes) are
stable chimeras, recorded for higher $\displaystyle R$ values. Other parameter
values producing interesting stripe chimeras in this area are the pairs
$\displaystyle(\varphi,\>\>R)$: $\displaystyle(300^{\circ},\>\>14)$,
$\displaystyle(299^{\circ},\>\>15))$, $\displaystyle(295^{\circ},\>\>16)$,
$\displaystyle(298^{\circ},\>\>19)$, $\displaystyle(299^{\circ},\>\>19)$,
$\displaystyle(296^{\circ},\>\>20)$, $\displaystyle(297^{\circ},\>\>20)$,
$\displaystyle(298^{\circ},\>\>20)$ and $\displaystyle(299^{\circ},\>\>20)$.

Figure 4:  Synchronous and asynchronous stripes for high
$\displaystyle\varphi$-values. Spatial $\displaystyle X-$profiles for: a)
$\displaystyle(\varphi=298^{\circ},\>\>R=15)$, b)
$\displaystyle(\varphi=299^{\circ},\>\>R=14)$, c)
$\displaystyle(\varphi=295^{\circ},\>R=20)$ and d)
$\displaystyle(\varphi=300^{\circ},\>R=19)$. Other parameters as in Fig. 2.
All simulations start from the same random initial conditions. A related video
is presented in the Supplementary Material (Multimedia view corresponding to
image Fig. 4c (movie_fig04c_brusselator_sig_1.00_phi_295_R_20.mp4)).

We can also observe transitions from coherent traveling linear fronts to
traveling chimera cores via the production of solitary states for rotation
angle values around $\displaystyle\varphi\approx 220^{\circ}-240^{\circ}$, in
the center of the parameter regions of Fig. 2. In this parameter region
$\displaystyle Z<1$ and $\displaystyle\Delta\omega>0$. Keeping $\displaystyle
R=21$, we gradually increase $\displaystyle\varphi$ and we present the
transitions from traveling coherent fronts to traveling coherent/incoherent
domains, in Fig. 4. In this figure, the spatial $\displaystyle X-$profiles for
$\displaystyle\varphi=220^{\circ},\>R=21$ in panel a show the development of
coherent traveling fronts. As the rotation coupling angle increases to
$\displaystyle\varphi=225^{\circ}$ we note in panel b the presence of a few
(three) solitary Brusselators whose state variables escape from those of the
neighbors. When $\displaystyle\varphi$ is further increased to
$\displaystyle\varphi=228^{\circ}$, more solitaries mobilize causing
perturbations in the linear structure of the traveling fronts, see panel c.
For $\displaystyle\varphi=229^{\circ}$ the solitaries organize into incoherent
regions as shown in panel d. Further increase of $\displaystyle\varphi$ to
$\displaystyle\varphi=236^{\circ}$ causes increase in the size of the
incoherent regions and total deformation of the linear structure of the
traveling fronts. Finally, for $\displaystyle\varphi=236^{\circ}$ the coherent
fronts break totally and incoherent cores are formed as satellites to
traveling coherent domains. Indeed, one can observe ramifications in the
$\displaystyle\Delta\omega$ values in the central part of Fig. 2b, which
justify the presence of traveling front (coherent structures) with
$\displaystyle\Delta\omega=0$ (black regions) next to traveling chimera states
with $\displaystyle\Delta\omega>0$ (red and blue regions). Similar behavior
involving transitions from coherent traveling fronts to moving chimera states
via solitary mobilization is observed for large and intermediate values of
$\displaystyle R$ and intermediate $\displaystyle\varphi$ domains.

Figure 5:  Transitions from coherent stripes to solitaries and further to
coherent and incoherent traveling domains. Spatial $\displaystyle X-$profiles
for: a) $\displaystyle(\varphi=220^{\circ},\>\>R=21)$, b)
$\displaystyle(\varphi=225^{\circ},\>\>R=21)$, c)
$\displaystyle(\varphi=228^{\circ},\>\>R=21)$, d)
$\displaystyle(\varphi=229^{\circ},\>\>R=21)$, e)
$\displaystyle(\varphi=230^{\circ},\>\>R=21)$ and f)
$\displaystyle(\varphi=236^{\circ},\>\>R=21)$. Other parameters as in Fig. 2.
All simulations start from the same random initial conditions. A related video
is presented in the Supplementary Material (Multimedia view corresponding to
image Fig. 5e (movie_fig05e_brusselator_sig_1.00_phi_230_R_21.mp4)).

In the vicinity of $\displaystyle\varphi\approx\pi(=180^{\circ})$ we can
observe multiple coherent spots surrounded by an incoherent halos. Depending
on the parameter value, the size and number of the coherent spots varies to
cover exhaustively the 2D lattice. Representative examples are provided in
Fig. 6. In panel a, for $\displaystyle\varphi=160^{\circ},\>\>R=16$, three
zones of interconnected coherent spots are observed. These zones are separated
by coherent ones of varying width. The borders between the coherent spots and
the coherent zones consist of incoherent elements. The coherent spots and
coherent zones are characterized by a small, but persistent phase difference.
When the connectivity range increases to $\displaystyle R=21$ in panel b, the
number (and form) of the coherent spots decreases (from 9 to 4) and the area
of the incoherent domains grows. The coherent strings also degenerate into 4
coherent domains lagging behind with respect to the value of the state
variable as compared to the coherent spots. In panel c, for different
parameter values (in the same parameter domain)
$\displaystyle\varphi=165^{\circ},\>\>R=19$, the size of the connected
coherent spots increases and their multiplicity decreases to 4. Here the phase
difference between the coherent spots and the coherent zones is considerable,
the former being in the red-yellow colors (highest $\displaystyle X$-values),
while the latters are in the black-blue colors (lowest $\displaystyle
X$-values). By increasing slightly $\displaystyle\varphi$ to $\displaystyle
166^{\circ}$, in panel d, the phase difference between the coherent spots and
the coherent zones vanishes and the connecting incoherent domains are now
visible in vivid colors.

Figure 6:  Multiple connected and disconnected coherent spots. Spatial 2D
$\displaystyle X-$profiles for: a)
$\displaystyle(\varphi=160^{\circ},\>\>R=16)$, b)
$\displaystyle(\varphi=160^{\circ},\>\>R=21)$, c)
$\displaystyle(\varphi=165^{\circ},\>\>R=19)$ and d)
$\displaystyle(\varphi=166^{\circ},\>\>R=19)$. Other parameters as in Fig. 2.
All simulations start from the same random initial conditions.

Along the right border between turbulent states and coherence, stable
incoherent spots may develop for appropriate values of the system parameters.
These can be considered as shrinking or degeneration of the stripes which are
observed in nearby parameter areas (see Fig. 4). In Fig. 7 a variety of
incoherent spots are recorded, for parameter values $\displaystyle
302\leq\varphi\leq 304$ and different $\displaystyle R$.

Figure 7:  Single incoherent spots. Spatial 2D $\displaystyle X-$profiles for:
a) $\displaystyle(\varphi=303^{\circ},\>\>R=6)$, b)
$\displaystyle(\varphi=302^{\circ},\>\>R=10)$, c)
$\displaystyle(\varphi=304^{\circ},\>\>R=8)$ and d)
$\displaystyle(\varphi=302^{\circ},\>\>R=18)$. Other parameters as in Fig. 2.
All simulations start from the same random initial conditions. Two related
videos are presented in the Supplementary Material (Multimedia view
corresponding to image Fig. 7b
(movie_fig07b_brusselator_sig_1.00_phi_302_R_10.mp4) and Fig. 7c
(movie_fig07c_brusselator_sig_1.00_phi_304_R_08.mp4)).

Besides the patterns discussed above, we also report the following ones:

  * •

For $\displaystyle\varphi\approx 170^{\circ}-200^{\circ}$, traveling square
coherent domains develop. Depending on the values of $\displaystyle R$ the
coherent squares can be surrounded by incoherent halos.

  * •

Linear coherent stripes cross the system for large values of $\displaystyle
R$, e.g., $\displaystyle R=21$ and $\displaystyle 203^{\circ}\leq\varphi\leq
225^{\circ}$ and other parameters values in the same area.

  * •

Combinations of coherent linear parts and incoherent irregular shapes are
often recorded. These are observed for smaller values of $\displaystyle R$,
e.g., $\displaystyle R=6$ and similar values.

  * •

Scroll waves are recorded for parameter values near the stable chimera
patterns, when the incoherent domains dissolve, e.g., for
$\displaystyle\varphi\approx 300^{\circ}$ and $\displaystyle 2\leq R\leq 10$.
Examples of these and other turbulent chimera patterns are shown in Appendix
B.

  * •

Results for the case $\displaystyle R=1$, which corresponds to the diffusion
limit, are also shown in Appendix B.

The parameter domain were the various patterns dominate were obtained for the
working parameter set given in the first paragraph of the present section. For
different generic parameter values the borders between the various domains of
coherence are adjusted appropriately.

##  V General comments on spatial and spatiotemporal pattern formation in the
Brusselator model

From the results in the previous two sections, it becomes clear that in
different parameter domains the Brusselator network exhibits different spatial
(Turing) and spatiotemporal (chimera) patterns of high complexity.

Indeed, Turing patterns describe in section III are spatial patterns which are
stable in time after a transient period, while in chimera states all elements
perform nonlinear oscillations, but they are spatially organized in coherent
and incoherent domains. As an additional evidence of the intrinsic difference
in the temporal evolution between these two types of patterns, we present in
the two figures below the temporal evolution of particular oscillators, in
Fig. 8, and their corresponding Fourier spectra, in Fig. 9. We have chosen to
demonstrate this for parameters as in Fig. 1a1 corresponding to the Turing
regime and Fig. 4c in the chimera regime.

Figure 8:  Temporal evolution of Brusselator elements: a) For parameter values
in the Turing regime, as in Fig. 1a1. The black solid (lower) line corresponds
to the element $\displaystyle(32,7)$ in the center of a spot and the red
(upper) line to the element $\displaystyle(32,14)$ outside the spot areas.
Other parameters as in Fig. 1a1. b) For parameter values in the chimera
regime, as in Fig. 4c. The red solid line corresponds to the element
$\displaystyle(32,60)$ in the coherent region and the black line to the
element $\displaystyle(32,32)$ in the incoherent stripe. Other parameters as
in Fig. 4c.

From Fig. 8a, it is evident that the Brusselators with parameters set at the
Turing regime reach a stable steady state, different for each element. Namely,
the $\displaystyle X-$variable of Brusselator at position (32,7), which
corresponds to the center of a spot in Fig. 1a1, reaches value $\displaystyle
X(32,7)\to\sim 3$, while the one at position (32,14), which is located outside
the spot areas, reaches $\displaystyle X(32,14)\to\sim 5$. To the contrary,
when the Brusselators are set in the chimera regime, they all oscillate in
time, as shown in Fig. 8b. In Fig. 8b (red line), we have chosen to
demonstrate the temporal evolution of a node in the coherent domain, (node
(32,60) in Fig. 4c). In the same figure the black line shows the evolution of
a node in the incoherent domain (node (32,32) in Fig. 4c). Both elements show
similar oscillations, however, we may observe that the Brusselator which
belongs to the incoherent region, position (32,32), is characterized by lower
frequency of oscillations than the one at position (32,60) in the coherent
region. This is a well known property of chimera states, where the coherent
and incoherent regions are characterized by different frequenciesKuramoto and
Battogtokh (2002); Omelchenko  _et al._ (2015); Schöll (2016); Omel’Chenko
(2018), as will be also verified in the next figure. Similar results and
conclusions are obtained for all other parameters in the chimera regime.

Figure 9:  Fourier spectra corresponding to the temporal evolution of the
Brusselator elements shown in Fig. 8b. The parameter values and line
descriptions follow the ones in Fig. 8b. For the extraction of the Fourier
spectra DFT Fourier transform was used and time series covered approximately
approximately 170 cycles, in both cases, excluding the first 50 cycles as
transient.

To investigate the difference in the frequencies in the coherent and
incoherent domains of the chimera states, as an example we plot in Fig. 9 the
Fourier spectra of the two time series shown in Fig. 8b. We choose again the
same two Brusselator elements at positions (32,60), in the coherent region,
plotted with the red line and at position (32,32), in the incoherent region,
plotted with black line. From the plot, it is evident that for both elements
only one frequency prevails, together with the corresponding harmonics. The
maxima in the two curves are at different positions, consistently for the main
frequencies and their harmonics. This observation corroborates the statement
in the previous paragraph on the difference in the oscillator frequencies in
the coherent and the incoherent regions of the chimera states.

##  VI Conclusions

We have investigated the conditions under which Turing patterns and chimera
states can be produced in the same nonlinear dynamical system. Using as
exemplary case the Brusselator model, we have demonstrated Turing patterns for
parameter regions where the single node dynamics presents fixed points. The
same dynamical system exhibits chimera states if the parameters of the single
element are chosen in the oscillatory (limit cycle) regime. Connectivity-wise,
for the emergence of the Turing patterns small connectivity neighborhoods are
required, mimicking the nearest neighbor, Laplacian diffusion, while for the
production of chimera states the connectivity regions need to be of
intermediate range.

Based on the above results and discussions, we may add that it is not possible
to obtain chimera states and Turing patterns for the same values of parameters
$\displaystyle A$ and $\displaystyle B$ in the Brusselator model, evenif the
other (connectivity) parameters are different. This is because for the
emergence of Turing patterns (stable steady states) the single element
dynamics need to lead to fixed points, while for the emergence of chimera
states (oscillatory steady states) a limit cycle is necessary in the single
element dynamics. The above conclusions hold provided that the coupling terms
do not drastically alter the nodal dynamics. Additional numerical and
analytical studies using different coupled nonlinear models are necessary to
further elucidate the necessary and sufficient conditions for the presence of
chimera states in these systems.

###### Acknowledgements.

The author (A.P.) acknowledges helpful discussions with K. Anesiadis. This
work was supported by computational time granted from the Greek Research &
Technology Network (GRNET) in the National HPC facility - ARIS, under project
IDs PR009012 and PR12015.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Discrete form of Brusselator with diffusion

In this appendix we discretize the Brusselator model with diffusion and bring
it to the general form, Eq. (4), in order to use a generic system of
differential equations, both for the production of Turing patterns and chimera
states (at different regions of the parameter space).

Starting from Eq. (10), we specify the location $\displaystyle(i,j)$ of a
Brusselator node and analyze the Laplacian operator around this location.
Then, Eq. (10) reduces to:

$\displaystyle\displaystyle\frac{dX_{ij}}{dt}$
$\displaystyle\displaystyle=A-(B+1)X_{ij}+X_{ij}^{2}Y_{ij}+D_{X}\nabla^{2}X_{ij}$
(13a) $\displaystyle\displaystyle\frac{dY_{ij}}{dt}$
$\displaystyle\displaystyle=BX_{ij}-X_{ij}^{2}Y_{ij}+D_{Y}\nabla^{2}Y_{ij}.$
(13b)

where $\displaystyle\nabla^{2}$ is the Laplacian operator in two dimensions
and $\displaystyle D_{X}$, $\displaystyle D_{Y}$ are the diffusion constants
of the $\displaystyle X$ and $\displaystyle Y$ variables, respectively. By
expanding the Laplacian operator in space around position $\displaystyle(i,j)$
we obtain:

$\displaystyle\displaystyle\frac{dX_{ij}}{dt}$
$\displaystyle\displaystyle=A-(B+1)X_{ij}+X_{ij}^{2}Y_{ij}+\frac{D_{X}}{(\Delta
x)^{2}}\left[X_{i+1,j}-2X_{ij}+X_{i-1,j}\right]+\frac{D_{X}}{(\Delta
y)^{2}}\left[X_{i,j+1}-2X_{ij}+X_{i,j-1}\right]$ (14a)
$\displaystyle\displaystyle\frac{dY_{ij}}{dt}$
$\displaystyle\displaystyle=BX_{ij}-X_{ij}^{2}Y_{ij}+\frac{D_{Y}}{(\Delta
x)^{2}}\left[Y_{i+1,j}-2Y_{ij}+Y_{i-1,j}\right]+\frac{D_{Y}}{(\Delta
y)^{2}}\left[Y_{i,j+1}-2Y_{ij}+Y_{i,j-1}\right].$ (14b)

In Eqs. (14), $\displaystyle\Delta x$ and $\displaystyle\Delta y$ represent
the infinitesimal displacement on the x- and y-axis, respectively, as per the
derivative definition. In the case of discretization, the smallest possible
displacement on the lattice is one lattice site, and therefore
$\displaystyle\Delta x=\Delta y=h=1$ in the network discretization
approximation. Seeing from a different point of view, the displacement values
$\displaystyle\Delta x$ and $\displaystyle\Delta y$ can be absorbed in the
factor $\displaystyle D_{X}$ and $\displaystyle D_{Y}$. Gathering together the
terms which are linked to the neighboring nodes, Eq. (14) reduce further to:

$\displaystyle\displaystyle\frac{dX_{ij}}{dt}$
$\displaystyle\displaystyle=A-(B+1)X_{ij}+X_{ij}^{2}Y_{ij}+\frac{D_{X}}{(\Delta
x)^{2}}\sum_{k=i-1}^{i+1}\left[X_{kj}-X_{ij}\right]+\frac{D_{X}}{(\Delta
y)^{2}}\sum_{l=j-1}^{j+1}\left[X_{il}-X_{ij}\right]$ (15a)
$\displaystyle\displaystyle\frac{dY_{ij}}{dt}$
$\displaystyle\displaystyle=BX_{ij}-X_{ij}^{2}Y_{ij}+\frac{D_{Y}}{(\Delta
x)^{2}}\sum_{k=i-1}^{i+1}\left[Y_{kj}-Y_{ij}\right]+\frac{D_{Y}}{(\Delta
y)^{2}}\sum_{l=j-1}^{j+1}\left[Y_{il}-Y_{ij}\right].$ (15b)

Finally, gathering the terms on the right-hand-sides and because
$\displaystyle\Delta x=\Delta y=h=1$, Eqs. (15) result in:

$\displaystyle\displaystyle\frac{dX_{ij}}{dt}$
$\displaystyle\displaystyle=A-(B+1)X_{ij}+X_{ij}^{2}Y_{ij}+\frac{D_{X}}{h^{2}}\sum_{\\{k,l\\}}\left[X_{kl}-X_{ij}\right]$
(16a) $\displaystyle\displaystyle\frac{dY_{ij}}{dt}$
$\displaystyle\displaystyle=BX_{ij}-X_{ij}^{2}Y_{ij}+\frac{D_{Y}}{h^{2}}\sum_{\\{k,l\\}}\left[Y_{kl}-Y_{ij}\right].$
(16b)

Equations (16) have the same form with the more general Eqs. (4), with the
following restrictions:

$\displaystyle\displaystyle c_{ij;kl}\>w_{XX}=D_{X}/h^{2},$ (17a)
$\displaystyle\displaystyle c_{ij;kl}\>w_{YY}=D_{Y}/h^{2},$ (17b)
$\displaystyle\displaystyle w_{XY}=w_{YX}=0,$ (17c)
$\displaystyle\displaystyle R=1.$ (17d)

I.e., the more generic form of the Brusselator network, Eqs. (4), can give
rise to Turing patterns when cross-coupling terms are absent from the dynamics
and in the limit of local, nearest-neighbor exchanges, i.e., under the
conditions Eqs. (17).

##  Appendix B Additional turbulent chimeric patterns

In the present section additional interesting patterns are recorded, composed
by coexisting coherent and incoherent traveling patterns.

First, in Fig. 10 turbulent traveling patterns are shown. Of particular
interest is panel d, which consists of traveling fronts (red colors). One of
the fronts forms a spiral pattern at the interior of which an incoherent core
dominates. Similar incoherent cores are reported in Refs. Totz  _et al._
(2018, 2020) where the ZBKE model is used.

Figure 10:  Turbulent patterns. Spatial 2D $\displaystyle X-$profiles for: a)
$\displaystyle(\varphi=303^{\circ},\>\>R=6)$, b)
$\displaystyle(\varphi=302^{\circ},\>\>R=10)$, c)
$\displaystyle(\varphi=304^{\circ},\>\>R=8)$ and d)
$\displaystyle(\varphi=302^{\circ},\>\>R=18)$. Other parameters as in Fig. 2.
All simulations start from the same random initial conditions.

Of special interest are the patterns developed for $\displaystyle R=1$ a nd
with various $\displaystyle\varphi$ values. Note that the case $\displaystyle
R=1$ does not strictly correspond to the discretized regular diffusion, due to
the influence of the rotational matrix $\displaystyle\varphi$. Only the case
$\displaystyle\varphi=0,\>\>R=1$ can be considered as regular diffusion. The
results for incremental values of $\displaystyle\varphi$ are shown in Fig.
11a-h.

Figure 11:  Spatiotemporal patterns for $\displaystyle R=1$. Spatial 2D
$\displaystyle X-$profiles for: a) $\displaystyle\varphi=73^{\circ}$, b)
$\displaystyle\varphi=74^{\circ}$, c) $\displaystyle\varphi=120^{\circ}$, d)
$\displaystyle\varphi=180^{\circ}$, e) $\displaystyle\varphi=200^{\circ}$, f)
$\displaystyle\varphi=310^{\circ}$, g) $\displaystyle\varphi=314^{\circ}$ and
h) $\displaystyle\varphi=315^{\circ}$. The x-axis correspond to the horizontal
arrangement of the Brusselators and the y-axis to their vertical arrangement
on the 2D square lattice network. Other parameters as in Fig. 2. All
simulations start from the same random initial conditions.

For small $\displaystyle\varphi\leq 73^{\circ}$ and large
$\displaystyle\varphi\geq 314^{\circ}$ values of the rotation angle, the
system presents coherence as shown in Fig. 11a and h. Between these values a
certain degree of incoherence reins in the system as presented in panels b to
g of the same figure. The borderline values written here were obtained for the
working parameter set given in the first paragraph of section IV. For
different generic parameter values the borders between full coherence and
partial coherence (chimera states) are appropriately adjusted.

Finally, a notable stable pattern was recorded in the vicinity of the
parameter regions where the coherent stripes were found, for $\displaystyle
87^{\circ}\leq\varphi\leq 92^{\circ}$ and $\displaystyle 5\leq R\leq 12$. An
example is provided in Fig. 12 for ($\displaystyle\varphi=91^{\circ},\>R=10$);
it consists of two asymmetric coherent spots. The asymmetric spots oscillate
out of phase, as one can see from the present snapshot where at the same time
the $\displaystyle X$-variable of the largest spot is $\displaystyle\approx
0.5$ (blue color) and $\displaystyle X\approx 3$ (orange color) for the
smaller one. Traces of similar patterns are also seen in the surroundings of
this parameter value.

Figure 12:  Spatial 2D $\displaystyle X-$profiles for the two asymmetric
stable coherent spots occurring at :
$\displaystyle(\varphi=91^{\circ},\>\>R=10)$. Other parameters as in Fig. 2.
All simulations start from the same random initial conditions.

## References

  * Epstein and Pojman (1998) I. Epstein and J. Pojman,  _An introduction to nonlinear chemical dynamics_ (Oxford University Press, Oxford, 1998).
  * Turing (1952) A. Turing, Philosophical Transactions of the Royal Society B 237, 37 (1952).
  * Nicolis and Prigogine (1977) G. Nicolis and I. Prigogine,  _Self-organization in nonequilibrium systems_ (Wiley, New York, 1977).
  * Lefever and Nicolis (1971) R. Lefever and G. Nicolis, Journal of Theoretical Biology 30, 267 (1971).
  * Nicolis (1995) G. Nicolis,  _Introduction to Nonlinear Science_ (Cambridge University Press, Cambridge, 1995).
  * Epstein and Showalter (1996) I. Epstein and K. Showalter, Journal of Physical Chemistry 100, 13132 (1996).
  * Mikhailov and Showalter (2006) A. Mikhailov and K. Showalter, Physics Reports 425, 79 (2006).
  * Winfree (1978) A. T. Winfree, Progress in Theoretical Chemistry 4, 1 (1978).
  * Kepper  _et al._ (1991) P. D. Kepper, V. Castets, E. Dulos, and J. Boissonade, Physica D 49, 161 (1991).
  * Castets  _et al._ (1990) V. Castets, E. Dulos, J. Boissonade, and P. D. Kepper, Physical Review Letters 64, 2953 (1990).
  * Borckmans  _et al._ (1987) P. Borckmans, G. Dewel, D. Walgraef, and Y. Katayama, Journal of Statistical Physics 48, 1031 (1987).
  * Verdasca  _et al._ (1992) J. Verdasca, A. D. Wit, G. Dewel, and P. Borckmans, Physics Letters A 168, 194 (1992).
  * Borckmans  _et al._ (1992) P. Borckmans, A. D. Wit, and G. Dewel, Physica A 188, 137 (1992).
  * Pena and Perez-Garcia (2001) B. Pena and C. Perez-Garcia, Physical Review E 64, 056213 (2001).
  * Field  _et al._ (1972) R. J. Field, E. Koros, and R. M. Noyes, Journal of the Americal Chemical Society 94, 8649 (1972).
  * Field and Noyes (1974) R. J. Field and R. M. Noyes, Journal of Chemical Physics 60, 1877 (1974).
  * Becker and Field (1985) P. K. Becker and R. J. Field, Journal of Chemical Physics 89, 118 (1985).
  * Berenstein  _et al._ (2014) I. Berenstein, D. Bullara, and Y. D. Decker, Chaos 24, 033129 (2014).
  * Hildebrand  _et al._ (2001) M. Hildebrand, A. Skødt, and K. Showalter, Physical Review Letters 87, 088303 (2001).
  * Zhabotinsky  _et al._ (1993) A. M. Zhabotinsky, F. Buchholtz, A. B. Kiyatkin, and I. R. Epstein, Journal of Physical Chemistry 97, 75 (1993).
  * Nkomo  _et al._ (2013) S. Nkomo, M. R. Tinsley, and K. Showalter, Physical Review Letters 110, 244102 (2013).
  * Kuramoto and Battogtokh (2002) Y. Kuramoto and D. Battogtokh, Nonlinear Phenomena in Complex Systems 5, 380 (2002).
  * Kuramoto (2002) Y. Kuramoto, in  _Nonlinear Dynamics and Chaos: Where do we go from here?_ , edited by S. J. Hogan, A. R. Champneys, A. R. Krauskopf, M. di Bernado, R. E. Wilson, H. M. Osinga, and M. E. Homer (CRC Press, 2002) pp. 209–227.
  * Abrams and Strogatz (2004) D. M. Abrams and S. H. Strogatz, Physical Review Letters 93, 174102 (2004).
  * Bolotov  _et al._ (2016) M. I. Bolotov, G. V. Osipov, and A. Pikovsky, Physical Review E 93, 032202 (2016).
  * Goldschmidt  _et al._ (2019) R. J. Goldschmidt, A. Pikovsky, and A. Politi, Chaos 29, 071101 (2019).
  * Bolotov  _et al._ (2020) M. I. Bolotov, L. A. Smirnov, G. V. Osipov, and A. Pikovsky, Physical Review E 102, 042218 (2020).
  * Sakaguchi (2006) H. Sakaguchi, Physical Review E 73, 031907 (2006).
  * Omelchenko  _et al._ (2013) I. Omelchenko, O. E. Omel’chenko, P. Hövel, and E. Schöll, Physical Review Letters 110, 224101 (2013).
  * Omelchenko  _et al._ (2015) I. Omelchenko, A. Provata, J. Hizanidis, E. Schöll, and P. Hövel, Physical Review E 91, 022917 (2015).
  * Schmidt  _et al._ (2017) A. Schmidt, T. Kasimatis, J. Hizanidis, A. Provata, and P. Hövel, Physical Review E 95, 032224 (2017).
  * Zakharova  _et al._ (2017) A. Zakharova, N. Semenova, V. S. Anishchenko, and E. Schöll, Chaos 27, 114320 (2017).
  * Shepelev and Vadivasova (2019) I. A. Shepelev and T. E. Vadivasova, Communications in Nonlinear Science and Numerical Simulation 79, 104925 (2019).
  * Hizanidis  _et al._ (2014) J. Hizanidis, V. Kanas, A. Bezerianos, and T. Bountis, International Journal of Bifurcations and Chaos 24, 1450030 (2014).
  * Ulonska  _et al._ (2016) S. Ulonska, I. Omelchenko, A. Zakharova, and E. Schöll, Chaos 26, 094825 (2016).
  * Yeldesbay  _et al._ (2014) A. Yeldesbay, A. Pikovsky, and M. Rosenblum, Physical Review Letters 112, 144103 (2014).
  * Gjurchinovski  _et al._ (2017) A. Gjurchinovski, E. Schöll, and A. Zakharova, Physical Review E 95, 042218 (2017).
  * Tumash  _et al._ (2018) L. Tumash, E. Panteley, A. Zakharova, and E. Schöll, European Physical Journal B 92, 100 (2018).
  * Olmi  _et al._ (2015) S. Olmi, E. A. Martens, S. Thutupalli, and A. Torcini, Physical Review E 92, 030901 (2015).
  * Tsigkri-DeSmedt  _et al._ (2017) N. D. Tsigkri-DeSmedt, J. Hizanidis, E. Schöll, P. Hövel, and A. Provata, The European Physical Journal B 90, 139 (2017).
  * Tsigkri-DeSmedt  _et al._ (2018) N. D. Tsigkri-DeSmedt, I. Koulierakis, G. Karakos, and A. Provata, The European Physical Journal B 91, 305 (2018).
  * Maistrenko  _et al._ (2017) V. Maistrenko, O. Sudakov, O. Osiv, and Y. Maistrenko, European Physical Journal: Special Topics 226, 1867 (2017).
  * Laing (2017) C. R. Laing, SIAM Journal on Applied Dynamical Systems 16, 974 (2017).
  * Omel’Chenko and Knobloch (2019) O. E. Omel’Chenko and E. Knobloch, New Journal of Physics 21, 093034 (2019).
  * Maistrenko  _et al._ (2015) Y. Maistrenko, O. Sudakov, O. Osiv, and V. Maistrenko, New Journal of Physics 17, 073037 (2015).
  * Kasimatis  _et al._ (2018) T. Kasimatis, J. Hizanidis, and A. Provata, Physical Review E 97, 052213 (2018).
  * Maistrenko  _et al._ (2020) V. Maistrenko, O. Sudakov, and O. Osiv, Chaos 30, 063113 (2020).
  * Koulierakis  _et al._ (2020) I. Koulierakis, D. A. Verganelakis, I. Omelchenko, A. Zakharova, E. Schöll, and A. Provata, Chaos 30, 113137 (2020).
  * M.Hagerstrom  _et al._ (2012) A. M.Hagerstrom, T. E. Murphy, R. Roy, P. Hövel, I. Omelchenko, and E. Schöll, Nature Physics 8, 658 (2012).
  * Gambuzza  _et al._ (2014) L. V. Gambuzza, A. Buscarino, S. Chessari, L. Fortuna, R. Meucci, and M. Frasca, Physical Review E 90, 032905 (2014).
  * Martens  _et al._ (2013) E. A. Martens, S. Thutupalli, A. Fourrière, and O. Hallatschek, Proceedings of the National Academy of Sciences 110, 10563 (2013).
  * Kapitaniak and Kurths (2014) T. Kapitaniak and J. Kurths, European Physical Journal: Special Topics 223, 609 (2014).
  * Dudkowski  _et al._ (2016) D. Dudkowski, J. Grabski, J. Wojewoda, P. Perlikowski, Y. Maistrenko, and T. Kapitaniak, Scientific Reports 6, 29833 (2016).
  * Manoj  _et al._ (2021) K. Manoj, S. A. Pawar, and R. I. Sujith, Physical Review E 103, 022207 (2021).
  * Cherry and Fenton (2008) E. M. Cherry and F. H. Fenton, New Journal of Physics 10, 125016 (2008).
  * Mormann  _et al._ (2000) F. Mormann, K. Lehnertz, P. David, and C. E. Elger, Physica D 144, 358 (2000).
  * Mormann  _et al._ (2003) F. Mormann, T. Kreuz, R. G. Andrzejak, P. David, K. Lehnertz, and C. E. Elger, Epilepsy Res. 53, 173 (2003).
  * Santos  _et al._ (2017) M. S. Santos, J. D. Szezech, F. S. Borges, K. C. Iarosz, I. L. Caldas, A. M. Batista, R. L. Viana, and J. Kurths, Chaos, Solitons & Fractals 101, 86 (2017).
  * Tinsley  _et al._ (2012) M. R. Tinsley, S. Nkomo, and K. Showalter, Nature Physics 8, 662 (2012).
  * Wickramasinghe and Kiss (2013) M. Wickramasinghe and I. Z. Kiss, PLoS ONE 8 (2013).
  * Schmidt  _et al._ (2014) L. Schmidt, K. Schönleber, K. Krischer, and V. García-Morales, Chaos 24, 013102 (2014).
  * Rattenborg  _et al._ (2000) N. C. Rattenborg, C. J. Amlaner, and S. L. Lima, Neuroscience and Biobehavioral Reviews 24, 817 (2000).
  * Rattenborg (2006) N. C. Rattenborg, Naturwissenschaften 93, 413 (2006).
  * Ramlow  _et al._ (2019) L. Ramlow, J. Sawicki, A. Zakharova, J. Hlinka, J. C. Claussen, and E. Schöll, Europhysics Letters 126, 50007 (2019).
  * Andrzejak  _et al._ (2016) R. G. Andrzejak, C. Rummel, F. Mormann, and K. Schindler, Scientific Reports 6, 23000 (2016).
  * Ruzzene  _et al._ (2019) G. Ruzzene, I. Omelchenko, E. Schöll, A. Zakharova, and R. G. Andrzejak, Chaos 29, 0511031 (2019).
  * Bansal  _et al._ (2019) K. Bansal, J. O. Garcia, S. H. Tompson, T. Verstynen, J. M. Vettel, and S. F. Muldoon, Science Advances 5, eaau853 (2019).
  * Panaggio and Abrams (2015) M. J. Panaggio and D. Abrams, Nonlinearity 28, R67 (2015).
  * Schöll (2016) E. Schöll, European Physical Journal – Special Topics 225, 891 (2016).
  * Yao and Zheng (2016) N. Yao and Z. Zheng, International Journal of Modern Physics B 30, 1630002 (2016).
  * Majhi  _et al._ (2019) S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, Physics of Life Reviews 28, 100 (2019).
  * Omel’Chenko (2018) O. E. Omel’Chenko, Nonlinearity 31, R121 (2018).
  * Zakharova (2020) A. Zakharova,  _Chimera Patterns in Networks: Interplay between Dynamics, Structure, Noise, and Delay_ (Understanding Complex Systems, Springer International Publishing, 2020).
  * Omel’Chenko (2022) O. E. Omel’Chenko, Journal of Nonlinear Science 32, 22 (2022).
  * Taylor  _et al._ (2015) A. F. Taylor, M. R. Tinsley, and K. Showalter, Physical Chemistry Chemical Physics 17, 20047 (2015).
  * Nkomo  _et al._ (2016) S. Nkomo, M. R. Tinsley, and K. Showalter, Chaos 26, 094826 (2016).
  * Totz  _et al._ (2018) J. F. Totz, J. Rode, M. R. Tinsley, K. Showalter, and H. Engel, Nature Physics 14, 282 (2018).
  * Totz  _et al._ (2020) J. F. Totz, M. R. Tinsley, H. Engel, and K. Showalter, Scientific Reports 10, 7821 (2020).
  * Hizanidis  _et al._ (2015) J. Hizanidis, E. Panagakou, I. Omelchenko, E. Schöll, P. Hövel, and A. Provata, Physical Review E 92, 012915 (2015).
  * Franović  _et al._ (2021) I. Franović, O. E. Omel’Chenko, and M. Wolfrum, Physical Review E 104, L052201 (2021).
  * Hizanidis  _et al._ (2020) J. Hizanidis, N. Lazarides, and G. P. Tsironis, Chaos 30, 013115 (2020).
  * Kokolnikov  _et al._ (2006) T. Kokolnikov, T. Erneux, and J. Wei, Physica D 214, 63 (2006).
  * Budroni and Wit (2016) M. A. Budroni and A. D. Wit, Physical Review E 93, 062207 (2016).
  * Bullara  _et al._ (2013) D. Bullara, Y. D. Decker, and R. Lefever, Physical Review E 87, 062023 (2013).
  * Tildi  _et al._ (2016) M. Tildi, Y. Gandica, G. Sonnino, E. Averlant, and K. Panajotov, Entropy 18, 64 (2016).
  * Shoji  _et al._ (2007) H. Shoji, K. Yamada, D. Ueyama, and T. Ohta, Physical Review E 75, 046212 (2007).
  * Omel’Chenko  _et al._ (2010) O. E. Omel’Chenko, M. Wolfrum, and Y. L. Maistrenko, Physical Review E 81, 065201(R) (2010).
  * Mikhaylenko  _et al._ (2019) M. Mikhaylenko, L. Ramlow, S. Jalan, and A. Zakharova, Chaos 29, 023122 (2019).
  * Bick  _et al._ (2016) C. Bick, P. Ashwin, and A. Rodrigues, Chaos 26, 094814 (2016).


print(A.storage()) # Storage elements: [ 0 1 2 3 4 5 6 7 ]

# Make it non-contiguous

A.permute_(0,2,1)

print(A.is_contiguous()) # Output: False

# Note that the storage is not changed

print(A.storage()) # Storage elements: [ 0 1 2 3 4 5 6 7 ]

# Make it contiguous, thus, the elements are permuted in memory

A.contiguous_()

print(A.is_contiguous()) # Output: True

# Note that the storage now is changed

print(A.storage()) # Storage elements: [ 0 2 1 3 4 6 5 7 ]

Listing 77: Reorder elements of Storage after permute.

##### OpenMP.

If OpenMP is available, Cytnx uses it for the parallel execution of code on
CPUs. If Cytnx is compiled without OpenMP, Device.Ncpus (see Sec. 10.1) is
always 1. If OpenMP is enabled, one can change the environment variable
OMP_NUM_THREADS before executing the program to set a restriction on how many
threads the program can use. For example, the following code example makes
OpenMP and the internal functions of Cytnx use $16$ threads where possible:

⬇

# In the command prompt, before executing a Cytnx program:

export OMP_NUM_THREADS=16

Listing 78: Set the number of threads used by OpenMP.

The speedup when using OpenMP and eight threads instead of only one thread can
be seen in the benchmark in Fig. 18.

## References

  * [1] R. Orús,  _A practical introduction to tensor networks: Matrix product states and projected entangled pair states_ ,  Annals of Physics 349, 117 (2014),  10.1016/j.aop.2014.06.013.
  * [2] U. Schollwöck,  _The density-matrix renormalization group in the age of matrix product states_ ,  Annals of Physics 326(1), 96 (2011),  10.1016/j.aop.2010.09.012,  January 2011 Special Issue.
  * [3] J. I. Cirac, D. Pérez-García, N. Schuch and F. Verstraete,  _Matrix product states and projected entangled pair states: Concepts, symmetries, theorems_ ,  Rev. Mod. Phys. 93, 045003 (2021),  10.1103/RevModPhys.93.045003.
  * [4] F. Verstraete, T. Nishino, U. Schollwöck, M. C. Bañuls, G. K. Chan and M. E. Stoudenmire,  _Density matrix renormalization group, 30 years on_ ,  Nature Reviews Physics 5(5), 273 (2023),  10.1038/s42254-023-00572-5.
  * [5] H. Fehske, R. Schneider and A. Weiße,  _Computational Many-Particle Physics_ ,  Lecture Notes in Physics. Springer Berlin, Heidelberg,  ISBN 978-3-540-74685-0,978-3-642-09414-9,978-3-540-74686-7,  10.1007/978-3-540-74686-7 (2008).
  * [6] R. Orús,  _Tensor networks for complex quantum systems_ ,  Nature Reviews Physics 1(9), 538 (2019),  10.1038/s42254-019-0086-7.
  * [7] H. Weimer, A. Kshetrimayum and R. Orús,  _Simulation methods for open quantum many-body systems_ ,  Rev. Mod. Phys. 93, 015008 (2021),  10.1103/RevModPhys.93.015008.
  * [8] I. P. McCulloch,  _From density-matrix renormalization group to matrix product states_ ,  Journal of Statistical Mechanics: Theory and Experiment 2007(10), P10014 (2007),  10.1088/1742-5468/2007/10/P10014.
  * [9] G. K.-L. Chan and S. Sharma,  _The density matrix renormalization group in quantum chemistry_ ,  Annual Review of Physical Chemistry 62(1), 465 (2011),  10.1146/annurev-physchem-032210-103338,  PMID: 21219144.
  * [10] G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li and S. R. White,  _Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms_ ,  arXiv e-prints arXiv:1605.02611 (2016),  10.48550/arXiv.1605.02611.
  * [11] K. H. Marti and M. Reiher,  _The density matrix renormalization group algorithm in quantum chemistry_ ,  Zeitschrift für Physikalische Chemie 224(3-4), 583 (2010),  10.1524/zpch.2010.6125.
  * [12] G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li and S. R. White,  _Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms_ ,  The Journal of Chemical Physics 145(1), 014102 (2016),  10.1063/1.4955108.
  * [13] A. Baiardi and M. Reiher,  _The density matrix renormalization group in chemistry and molecular physics: Recent developments and new challenges_ ,  The Journal of Chemical Physics 152(4), 040903 (2020),  10.1063/1.5129672.
  * [14] L. Freitag and M. Reiher,  _The Density Matrix Renormalization Group for Strong Correlation in Ground and Excited States_ , chap. 7, pp. 205–245,  John Wiley & Sons, Ltd,  ISBN 9781119417774,  10.1002/9781119417774.ch7 (2020).
  * [15] Y. Cheng, Z. Xie and H. Ma,  _Post-density matrix renormalization group methods for describing dynamic electron correlation with large active spaces_ ,  The Journal of Physical Chemistry Letters 13(3), 904 (2022),  10.1021/acs.jpclett.1c04078.
  * [16] M. C. Bañuls and K. Cichy,  _Review on novel methods for lattice gauge theories_ ,  Reports on Progress in Physics 83(2), 024401 (2020),  10.1088/1361-6633/ab6311.
  * [17] Y. Meurice, R. Sakai and J. Unmuth-Yockey,  _Tensor lattice field theory for renormalization and quantum computing_ ,  Rev. Mod. Phys. 94, 025005 (2022),  10.1103/RevModPhys.94.025005.
  * [18] M. Dalmonte and S. Montangero,  _Lattice gauge theory simulations in the quantum information era_ ,  Contemporary Physics 57(3), 388 (2016),  10.1080/00107514.2016.1151199.
  * [19] M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik _et al._ ,  _Simulating lattice gauge theories within quantum technologies_ ,  The European Physical Journal D 74(8), 165 (2020),  10.1140/epjd/e2020-100571-8.
  * [20] E. Stoudenmire and D. J. Schwab,  _Supervised learning with tensor networks_ ,  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon and R. Garnett, eds., _Advances in Neural Information Processing Systems_ , vol. 29. Curran Associates, Inc.,  https://proceedings.neurips.cc/paper_files/paper/2016/file/5314b9674c86e3f9d1ba25ef9bb32895-Paper.pdf (2016).
  * [21] J. Liu, S. Li, J. Zhang and P. Zhang,  _Tensor networks for unsupervised machine learning_ ,  Phys. Rev. E 107, L012103 (2023),  10.1103/PhysRevE.107.L012103.
  * [22] S. Y.-C. Chen, C.-M. Huang, C.-W. Hsing and Y.-J. Kao,  _An end-to-end trainable hybrid classical-quantum classifier_ ,  Machine Learning: Science and Technology 2(4), 045021 (2021),  10.1088/2632-2153/ac104d.
  * [23] S. Y.-C. Chen, C.-M. Huang, C.-W. Hsing, H.-S. Goan and Y.-J. Kao,  _Variational quantum reinforcement learning via evolutionary optimization_ ,  Machine Learning: Science and Technology 3(1), 015025 (2022),  10.1088/2632-2153/ac4559.
  * [24] M. Kiffner and D. Jaksch,  _Tensor network reduced order models for wall-bounded flows_ ,  arXiv e-prints arXiv:2303.03010 (2023),  10.48550/arXiv.2303.03010.
  * [25] R. D. Peddinti, S. Pisoni, A. Marini, P. Lott, H. Argentieri, E. Tiunov and L. Aolita,  _Complete quantum-inspired framework for computational fluid dynamics_ ,  arXiv e-prints arXiv:2308.12972 (2023),  10.48550/arxiv.2308.12972.
  * [26] H.-J. Liao, K. Wang, Z.-S. Zhou, P. Zhang and T. Xiang,  _Simulation of ibm’s kicked ising experiment with projected entangled pair operator_ ,  arXiv e-prints arXiv:2308.03082 (2023),  10.48550/arXiv.2308.03082.
  * [27] T. Begušić and G. K.-L. Chan,  _Fast classical simulation of evidence for the utility of quantum computing before fault tolerance_ ,  arXiv e-prints arXiv:2306.16372 (2023),  10.48550/arXiv.2306.16372.
  * [28] J. Tindall, M. Fishman, M. Stoudenmire and D. Sels,  _Efficient tensor network simulation of ibm’s eagle kicked ising experiment_ ,  arXiv e-prints arXiv:2306.14887 (2023),  10.48550/arXiv.2306.14887.
  * [29] _Cytnx, Github repository_ ,  https://github.com/kaihsin/Cytnx.
  * [30] _Cytnx, documentation and user guide_ ,  https://kaihsinwu.gitlab.io/Cytnx_doc/Intro.html.
  * [31] M. Fishman, S. R. White and E. M. Stoudenmire,  _The ITensor Software Library for Tensor Network Calculations_ ,  SciPost Phys. Codebases p. 4 (2022),  10.21468/SciPostPhysCodeb.4.
  * [32] Y. Motoyama, T. Okubo, K. Yoshimi, S. Morita, T. Kato and N. Kawashima,  _Tenes: Tensor network solver for quantum lattice systems_ ,  Computer Physics Communications 279, 108437 (2022),  10.1016/j.cpc.2022.108437.
  * [33] J. Hauschild and F. Pollmann,  _Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy)_ ,  SciPost Phys. Lect. Notes p. 5 (2018),  10.21468/SciPostPhysLectNotes.5.
  * [34] _SyTen, tensor library_ ,  https://syten.eu/docs/.
  * [35] C. Roberts, A. Milsted, M. Ganahl, A. Zalcman, B. Fontaine, Y. Zou, J. Hidary, G. Vidal and S. Leichenauer,  _Tensornetwork: A library for physics and machine learning_ ,  arXiv e-prints arXiv:1905.01330 (2019),  10.48550/arXiv.1905.01330.
  * [36] J. Haegeman,  _Tensorkit.jl_ ,  https://jutho.github.io/TensorKit.jl.
  * [37] J. Gray,  _quimb: a python library for quantum information and many-body calculations_ ,  Journal of Open Source Software 3(29), 819 (2018),  10.21105/joss.00819.
  * [38] _Cytnx, API documentation_ ,  https://kaihsinwu.gitlab.io/cytnx_api/.
  * [39] _Cytnx, documentation and user guide, Section Linear algebra_ ,  https://kaihsinwu.gitlab.io/Cytnx_doc/guide/linalg.html.
  * [40] B. Pirvu, V. Murg, J. I. Cirac and F. Verstraete,  _Matrix product operator representations_ ,  New Journal of Physics 12(2), 025012 (2010),  10.1088/1367-2630/12/2/025012.
  * [41] S. Singh, R. N. C. Pfeifer and G. Vidal,  _Tensor network states and algorithms in the presence of a global u(1) symmetry_ ,  Phys. Rev. B 83, 115125 (2011),  10.1103/PhysRevB.83.115125.
  * [42] S. Singh, R. N. C. Pfeifer and G. Vidal,  _Tensor network decompositions in the presence of a global symmetry_ ,  Phys. Rev. A 82, 050301 (2010),  10.1103/PhysRevA.82.050301.
  * [43] S. Singh and G. Vidal,  _Tensor network states and algorithms in the presence of a global su(2) symmetry_ ,  Phys. Rev. B 86, 195114 (2012),  10.1103/PhysRevB.86.195114.
  * [44] T. Nishino and K. Okunishi,  _Corner transfer matrix renormalization group method_ ,  Journal of the Physical Society of Japan 65(4), 891 (1996),  10.1143/JPSJ.65.891.
  * [45] M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman and F. Verstraete,  _Faster methods for contracting infinite two-dimensional tensor networks_ ,  Phys. Rev. B 98, 235148 (2018),  10.1103/PhysRevB.98.235148.
  * [46] F. Verstraete and J. I. Cirac,  _Renormalization algorithms for quantum-many body systems in two and higher dimensions_ ,  arXiv e-prints cond-mat/0407066 (2004),  10.48550/arXiv.cond-mat/0407066.
  * [47] M. Schneider, J. Ostmeyer, K. Jansen, T. Luu and C. Urbach,  _Simulating both parity sectors of the hubbard model with tensor networks_ ,  Phys. Rev. B 104, 155118 (2021),  10.1103/PhysRevB.104.155118.
  * [48] M. Schneider,  _The Hubbard model on a honeycomb lattice with fermionic tensor networks_ ,  Ph.D. thesis, Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät,  10.18452/25393 (2022).
  * [49] _Cytnx, documentation and user guide, Section Iterative Solver_ ,  https://kaihsinwu.gitlab.io/Cytnx_doc/guide/itersol.html.
  * [50] S. Tomov, J. Dongarra and M. Baboulin,  _Towards dense linear algebra for hybrid GPU accelerated manycore systems_ ,  Parallel Computing 36(5-6), 232 (2010),  10.1016/j.parco.2009.12.005.
  * [51] _cuBLAS – Basic Linear Algebra on NVIDIA GPUs_ ,  https://developer.nvidia.com/cublas.
  * [52] _cuSOLVER – Direct Linear Solvers on NVIDIA GPUs_ ,  https://developer.nvidia.com/cusolver.
  * [53] _cuTENSOR – Tensor Linear Algebra on NVIDIA GPUs_ ,  https://developer.nvidia.com/cutensor.
  * [54] _cuQuantum – Accelerate quantum computing research._ ,  https://developer.nvidia.com/cublas.
  * [55] S. R. White,  _Density matrix formulation for quantum renormalization groups_ ,  Phys. Rev. Lett. 69, 2863 (1992),  10.1103/PhysRevLett.69.2863.
  * [56] G. Evenbly,  _Tensortrace: an application to contract tensor networks_ ,  10.48550/arXiv.1911.02558 (2019).
  * [57] H.-J. Liao, J.-G. Liu, L. Wang and T. Xiang,  _Differentiable programming tensor networks_ ,  Phys. Rev. X 9, 031041 (2019),  10.1103/PhysRevX.9.031041.
  * [58] J. Hasik, D. Poilblanc and F. Becca,  _Investigation of the Néel phase of the frustrated Heisenberg antiferromagnet by differentiable symmetric tensor networks_ ,  SciPost Phys. 10, 012 (2021),  10.21468/SciPostPhys.10.1.012.
  * [59] J. Naumann, E. Lennart Weerda, M. Rizzi, J. Eisert and P. Schmoll,  _variPEPS – a versatile tensor network library for variational ground state simulations in two spatial dimensions_ ,  arXiv e-prints arXiv:2308.12358 (2023),  10.48550/arXiv.2308.12358.
  * [60] A. Francuz, N. Schuch and B. Vanhecke,  _Stable and efficient differentiation of tensor network algorithms_ ,  arXiv e-prints arXiv:2311.11894 (2023),  10.48550/arXiv.2311.11894.
  * [61] P. Corboz, R. Orús, B. Bauer and G. Vidal,  _Simulation of strongly correlated fermions in two spatial dimensions with fermionic projected entangled-pair states_ ,  Phys. Rev. B 81, 165104 (2010),  10.1103/PhysRevB.81.165104.
  * [62] P. Corboz and G. Vidal,  _Fermionic multiscale entanglement renormalization ansatz_ ,  Phys. Rev. B 80, 165129 (2009),  10.1103/PhysRevB.80.165129.
  * [63] Z.-C. Gu,  _Efficient simulation of grassmann tensor product states_ ,  Phys. Rev. B 88, 115139 (2013),  10.1103/PhysRevB.88.115139.
  * [64] A. Yosprakob,  _GrassmannTN: A Python package for Grassmann tensor network computations_ ,  SciPost Phys. Codebases p. 20 (2023),  10.21468/SciPostPhysCodeb.20.
  * [65] F. Chollet _et al._ ,  _Keras_ ,  https://github.com/fchollet/keras (2015).