#### Stability mosaics in a forced Brusselator

THE EUROPEAN
PHYSICAL JOURNAL
Stability mosaics in a forced Brusselator
Joana G. Freire 1 2
Marcia R. Gallas 0 1 2 4
Jason A.C. Gallas 0 1 2 3 4
0 Complexity Sciences Center , 9225 Collins Avenue, Suite 1208, Surfside, FL 33154-3046 , USA
1 Departamento de F ́ısica, Universidade Federal da Para ́ıba , 58051-970 Joa ̃o Pessoa , Brazil
2 Instituto de Altos Estudos da Para ́ıba , Rua Silvino Lopes 419-2502, 58039-190 Jo ̃ao Pessoa , Brazil
3 Max-Planck-Institut fu ̈r Physik komplexer Systeme , 01187 Dresden , Germany
4 Institute for Multiscale Simulations, Friedrich-Alexander-Universita ̈t Erlangen-Nu ̈rnberg , 91052 Erlangen , Germany
We report a startling mosaic-like organization of stability phases found in the low-frequency limit of a driven Brusselator. Such phases correspond to periodic oscillations having a constant number of spikes per period. The mosaic is free from chaotic oscillations and is formed by an apparently infinite cascade of oscillations whose number of spikes grow without bound. Wide windows free from chaos but supporting unbounded quantities of complex oscillations are potentially of interest to operate driven oscillators such as lasers, electronic circuits, and biochemical pacemakers.
1 Introduction
In 1992, Lorenz remarked that “Now that computers have become ubiquitous,
carefully conceived numerical experiments can enable us to explore a fascinating
mathematical world that has not yet opened its doors to classical analytical
procedures” [1]. Nonlinear ordinary differential equations are ubiquitous in models of
lasers of various kinds, chemical and biochemical oscillators, cellular rhythms,
electronic circuits, mechanical oscillators, and several other natural phenomena [2, 3].
Studying the systematics of their dependence on control parameters arguably
qualifies as a fascinating mathematical world to explore and this is what we report in this
paper, for a periodically driven Brusselator [4].
Discriminating the nature of all possible oscillations and their stability is a
problem studied at least since the nineteenth century [5, 6] but still having many
unanswered questions [2, 3, 10]. At first, the work was essentially focused on studying
the phase-space, i.e., the space of variables of dynamical systems, with the help of
a number of approximate analytical methods. In the 1960s, numerical investigations
of the phase-space became more common and allowed one to go beyond the
limits of analytical approximations. Such investigations were restricted to a limited set
of parameters because the classification of the dynamical behaviors for a given set
presupposes, of course, the analysis of the full phase-space for it. Thus, numerical
investigation of large parameter domains is a rather computer demanding task.
Profiting from modern computer clusters, the purpose of this paper is to report
regularities discovered in stability diagrams of a periodically forced Brusselator, in the
low-frequency limit of the drive. A startling finding is that the low-frequency limit,
while characterized by wide regions completely free of chaos, contains systematic
accumulations of periodic oscillations whose waveforms have a number of spikes that
seem to grow without bound as a function of the control parameters.
2 The driven Brusselator
The system considered here is a periodically forced abstract chemical reaction known
as the Brusselator. The undriven version of the Brusselator was introduced in 1956 in
pioneering works dedicated to study undamped oscillations far from thermodynamic
equilibrium in open chemical systems governed by nonlinear kinetic laws [4]. The
driven Brusselator is defined by the equations [11–16]
= A − (B + 1)X + Y X2 + a cos(ωt),
= BX − Y X2.
Here, A and B are parameters controlling the chemical reaction, while a and ω are
two freely tunable parameters controlling the external drive whose impact we wish to
explore here.
Tomita, Kai, and Hikami [11] found two coexisting stable oscillations when
(A, B, a, ω) = (0.4, 1.2, 0.0018, 0.34). Subsequently, in a series of papers published
in 1978–82, Tomita and Kai reported additional results concerning the dynamics
observed when varying the frequency ω and amplitude a of the drive, when
maintaining the other two parameters fixed: A = 0.4 and B = 1.2. As typical in those
days, they were interested in the structure of resonance horns, in period-doubling
cascades ending in chaos, and in identifying quasi-periodic oscillations. Their
classical results are reviewed in reference [13]. The driven Brusselator was also studied by
Hao, Zhang, and coworkers in a series of papers published in 1982–89. They reported
detailed phase diagrams determined with the help of symbolic dynamics. Their
contributions are reviewed in references [14–16]. A recent review containing additional
stability diagrams is reference [4].
Although reporting a reasonable agreement between numerical and theoretical
results, Tomita, Kai, and Hikami remarked that despite the overall similarity between
computed plots and theoretical estimations, “ ... an appreciable discrepancy is noted
towards the low-frequency limit of entrainment.” As far as we known, the complicated
sequences of periodic oscillations unfolding in the low-frequency limit still remain to
be investigated. The low-frequency limit is of interest for applications, e.g., to control
the onset of pulsations, regular or not, in CO2 lasers with modulated losses [17–19] and
in semiconductor lasers [20, 21]. This is the oscillatory limit that we address here for
the Brusselator, showing that it harbors a plethora of unanticipated and remarkable
dynamical phenomena. In such limit the system is dominated by periodic oscillations
being free from chaotic motions. As detailed below, we find the global organization of
ω 0.1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Fig. 1. Absence of chaos in the low-frequency domain illustrated in two complementary
ways: (a) standard Lyapunov stability diagram, showing no λ > 0, and (b) isospike diagram
[17, 22–24] obtained by counting the number of spikes per period of X and representing them
in colors, as indicated. Lyapunov exponents are unable to discriminate individual oscillatory
phases while the isospike diagram makes their unfolding easily visible. The white box in
(b) reveals details of the wide-ranging mosaic structure, shown magnified in Figure 2. Each
panel displays the analysis of 600 × 600 parameter points. Here (a, B) = (0.4, 1.2).
the complex oscillations to deviate in several aspects from what is currently known.
Thus far, chaotic oscillations have dominated research in nonlinear oscillators.
However, after the coming-of-age of chaos it is becoming increasingly evident that the
global organization of periodic oscillations presents many unsuspected surprises and
that much still awaits to be discovered for them. No stability chart can be considered
complete if it does not register periodic and chaotic oscillations. In chemistry, periodic
oscillations already have an almost two-centuries long history, dating back at least
to the beginnings of the 19th century as evidenced by the many experimental works
collected and described in almost forgotten books [5, 6] and papers [7–9], a history
that deserves to be better known.
Before proceeding, we observe that the Brusselator is a “non-chemical” model
because of the termolecular step. However, despite this fact, it is widely used because
it gives the mass-action chemical expression of a very simple and fruitful model.
Furthermore, note that the periodic drive allows negative values of the variables. As
usual, following the literature, we are concerned with the dynamical aspects of the
model which, of course, describes just an abstract chemical reaction.
3 On representations of stability
Figure 1 illustrates the typical situation studied in this paper. For a given window in
the frequency versus the parameter A control space, it shows two distinct stability
diagrams: a typical Lyapunov stability diagram is shown in Figure 1a, while an isospike
diagram [17, 22–24] is shown in Figure 1b. This latter diagram represents in colors
parameter regions characterized by periodic oscillations having the same number of
spikes per period, as recorded in the colorbar under the diagram. Thus, the colors
allow one to visually grasp how the number of spikes varies when control parameters
are tuned. In the limit ω → 0, the isospike diagram shows a conspicuous symmetrical
accumulation of oscillations with an ever-increasing number of spikes, something not
visible in the Lyapunov diagram.
Both diagrams in Figure 1 were obtained by dividing the parameter window into
a mesh of N × N equally spaced parameter points and computing three quantities
at each point: the Lyapunov spectrum, the period of the oscillations (if any), and
the number of spikes contained in one period of every component of the periodic
oscillations. To this end, equation (1) was integrated numerically using the
standard fourth-order Runge-Kutta algorithm with fixed time-step h = 0.005. Integrations
were always started from an arbitrarily chosen initial condition, (X, Y ) = (0.16, 1.16),
that was maintained fixed for all integrations. The first 4 × 105 integration steps
were disregarded as transient-time needed to approach the attractor. The subsequent
80 × 105 time steps were used to determine the Lyapunov spectrum, Figure 1a.
After computing the Lyapunov spectrum, integrations were continued for an additional
40 × 105 time-steps, when we recorded up to 800 extrema (maxima and minima) of
the two variables of interest together with the instant when they occur. From this
record, we determined repetitions of the maxima, obtaining the period, if any, for
each parameter point.
As indicated by the colorbar in Figure 1b (and in similar Figures below), a palette
of 17 colors is used to represent the number of peaks (local maxima) per period.
Patterns having more than 17 peaks are plotted by periodically recycling the 17
basic colors “modulo 17”: Multiples 17 × of 17 are plotted with the color index
corresponding to (18 − ) mod 17, namely 34 ≡ 16 for = 2, 51 ≡ 15 for = 3, 68 ≡ 14
for = 4, etc. Note that substituting 17 by another number would only change the
order of colors in the diagram, not the size or boundaries of the phases that it contains.
In isospike diagrams, a black dot is normally plotted when no periodicity is
detected numerically. However, in Figure 1b, black represents a region where the
fixed integration step of the Runge-Kutta integrator becomes too big to correctly
detect all spikes in the oscillations. By decreasing the integration step one can eliminate
part of the black stripe until the step becomes too big again, needing to be reduced
once again, and so on. Such reductions make computations very time-consuming and
seem to add nothing essential to the diagram. In other words, the number of spikes
seems to increase continuously without bound when ω → 0.
The computation of stability diagrams is a numerically demanding task that
was performed using ad-hoc developed FORTRAN software to generate each figure
directly as Postscript bitmap output. Such computations have been described in
detail previously, e.g., in references [17, 25] where efficient methods to deal both with
numerical and experimental data are given. See also references [19, 26–28]. Isospike
diagrams require considerably less computations than Lyapunov diagrams. Of course,
they provide useful complementary information. Below, we focus on the diagrams that
are more adequate to our present goal.
4 Regular tiling of the low frequency limit
Figure 2 illustrates the regular tiling resulting from the rich variety of multipeaked
periodic oscillations in the ω → 0 limit, located in the region delimited by the white
box in Figure 1. In Figure 2, panels on the left column show diagrams obtained
when counting peaks of X, while panels on the right column display similar diagrams
but counting peaks of Y . Apart from a shift easily visible on the right portion of
Figures 2a–2b, both variables produce essentially the same “mosaic”. There,
periodic oscillations acquire extra spikes in a systematic way as parameters are changed
continuously.
An important feature in Figure 2 is the absence of chaotic phases. When moving
to the left on the panels, the number of spikes per period grows steadily by one. The
Fig. 2. Mosaic structures in the low-frequency domain observed by counting spikes per
period of the oscillations in X (left column) and Y (right column). White boxes in (a) and
(b) are shown magnified in (c) and (d), respectively. Resolution 600 × 600 parameter points.
Here (a, B) = (0.4, 1.2).
rapid color alternation indicates that the number of spikes per period grows very fast
towards ω = 0. As mentioned before, rather than chaos, the black region near ω = 0
reflects difficulties of the algorithm to separate series of very small amplitude peaks.
By reducing the integrator step and using longer transients we are able to complement
the diagram more and more near ω = 0. This makes us believe that, in fact, there is
no upper bound on the number of peaks, just on the amount of time and effort that
one is willing to spend in order to complete the diagram close to ω = 0. As may be
recognized from the figure, the double limit ω, A → 0 is even trickier and deserves
further study.
The upper panel in Figure 2 contains a number of representative points, indicated
by dots and numerical labels mn in the figure. Their coordinates, period and number
of peaks are listed in Table 1. The table also includes coordinates of points located
symmetrical with respect to the A = 0 axis.
Figures 3 and 4 depict the time-evolution of X and Y , respectively, for the
parameters listed in Table 1 and same initial conditions (X, Y ) = (0.16, 1.16) as before.
These figures show how the waveforms change when moving between the distinct
phases composing the mosaic. The sequence of waveforms reveals unexpected and
rather symmetrical evolutions. The phases containing the A = 0 axis correspond to
antiperiodic oscillations [29, 30], namely to oscillations obeying
s(t + T /2) = −s(t),
where the signal s(t) may be either s = X or s = Y , and T is the period of the
oscillation.
As may be appreciated from Figures 3 and 4, a more subtle symmetry emerges
between phases located symmetrically with respect to the A = 0 axis. Such phases
Fig. 4. Same
oscillations.
obey the relations:
Fig. 3. Systematic complexification of the X oscillatory patterns resulting in the symmetric
mosaic seen on the leftmost panel on the top. Here (a, B) = (0.4, 1.2). Oscillations 11, 32,
53 are antiperiodic (see text).
Figure 3, but showing the systematic complexification
s21(t + T /2) =
s31(t + T /2) =
s41(t + T /2) =
−s22(t),
−s33(t),
−s44(t),
peaks
1
2
2
3
3
3
4
4
4
4
5
5
5
5
5
peaks
1
2
2
3
3
3
4
4
4
4
5
5
5
5
5
s42(t + T /2) = −s43(t),
s51(t + T /2) = −s55(t),
s52(t + T /2) = −s54(t).
This shows antiperiodic oscillations s(t) to come in two distinct flavors: either as the
standard self-antiperiodic oscillations [29, 30], or as signals antiperiodic to another
oscillation in the group. As far as we know, this latter symmetry was not predicted
or observed before. The symmetry invariance of the equations of motion of the driven
Brusselator expose them, as well as the nice shift-symmetry underlying the
trigonometric drive. The mosaic reflects the symmetry (X, Y, A) → (−X, −Y, −A) and the
fact that cos(ωt + π) = − cos(ωt), underlying the T /2 shift.
Figures 3 and 4 allow one to understand the nature of the spike-adding mechanism
responsible for the regular tiling observed at low frequencies. Oscillations are
characterized by two distinct components: a slowly varying component together with fast
oscillations responsible for the sequence of small peaks. In the low-frequency limit,
the relatively slowly varying oscillations acquire localized trains of spikes that get
more and more low-amplitude spikes.
As mentioned before, while it is numerically demanding to repeatedly tune the
integrator and initial conditions to detect ever increasing number of spikes having tiny
differences in their amplitudes, numerical experiments show the scenario described to
remain consistent as the frequency gets smaller and smaller. In other words, we have
no reason to expect departures from the scenario described here for the low-frequency
limit.
What happens to the ω × A mosaic when the amplitude a of the external forcing
varies? This may be appreciated in Figure 5 for selected values of the amplitude, when
counting spikes of X. Counting spikes of Y results in similar diagrams, but with an
additional change of phases occurring along the line A = 0.425, and a few minor
details whose discussion is left for a future work. Essentially, the effect of increasing
a is to induce changes on the rightmost regions of the stability diagrams, basically
a strong compression, leaving the mosaic relatively the same, topologically invariant.
Fig. 5. Evolution of the stability mosaic as the amplitude a of the external drive
increases. Here, B = 1.2. Individual panels display the number of spikes per period of X
for (a) a = 0.2, (b) a = 0.3, (c) a = 0.35, (d) a = 0.4, (e) a = 0.45, (f) a = 0.5, (g) a = 0.7
and (h) a = 0.9. The structural organization of the leftmost region of the diagram remains
relatively unchanged while a marked evolution happens on rightmost region of the diagrams.
Note changes in scales. Resolution 600 × 600 parameter points.
A detailed characterization and explanation of the mosaics as a function of all control
parameters is an interesting but time consuming task.
5 Conclusions and outlook
This paper studied the behavior of the Brusselator driven by a low-frequency
periodic trigonometric force. This system was already studied with the basic aim
of understanding entrainment and the emergence of chaotic oscillations. Here, on
the contrary, we consider an interesting and broad complementary situation
concerning the global organization of stability phases resulting from periodic motions,
with an arbitrarily large number of spikes per period, in a region where chaos is
totally absent, namely the low-frequency limit. Of interest was to understand the
partitioning of large portions of the control parameter space into cascades of
stability phases and their intricate accumulations. The stability diagrams reported
here differ considerably from the ones familiar for, e.g., lasers [20, 21, 31],
chemical [27, 32], biochemical oscillators [33], and even for a novel psychological stress
variation model [34]. A startling observation reported here is an apparently
infinite sequence of signals that are not self-antiperiodic as usual but, instead, are
antiperiodic to another oscillation in the group. Another characteristic found in
the low-frequency limit is the presence of oscillations displaying changes in both
slow and fast scales. An interesting open question now is the investigation of the
low-frequency limit for other driven oscillators. The slow-fast scales of the
signals and their antiperiodic nature should not be difficult to observe in laboratory
experiments.
This work was supported by the Cluster of Excellence Engineering of Advanced
Materials, ZISC and FPS funded by the German Research Foundation (DFG) and by the
MaxPlanck Institute for the Physics of Complex Systems, Dresden, in the framework of the
Advanced Study Group on Optical Rare Events. JGF was supported by FCT, Portugal,
through the Post-Doctoral grant SFRH/BPD/101760/2014. JACG was supported by CNPq,
Brazil. The authors acknowledge use of the CESUP-UFRGS clusters located in Porto Alegre,
Brazil. Open access funding provided by Max Planck Society (Institut fu¨r Physik komplexer
Systeme).
1. E.N. Lorenz , J. Atmos. Sci. 49 , 2449 ( 1992 )
2. J. Argyris , G. Faust , M. Haase , R. Friedrich , An Exploration of Dynamical Systems and Chaos , 2nd edn. (Springer, New York, 2015 )
3. S. Strogatz , Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry, and Engineering, 2nd Edition, (Westview Press, Boulder, 2015 )
4. J.A.C. Gallas, Mod. Phys. Lett . B 29 , 1530018 ( 2015 ). Chapter in Peregrinations from Physics to Phylogeny: Essays on the Occasion of Hao Bailin's 80 Birthday , Edited by K.K. Phua , M. Ge (World Scientific, Singapore, 2016 )
5. R. Kremann , Die periodischen Erscheinungen in der Chemie (F . Enke Verlag, Stuttgart, 1913 )
6. E.S. Hedges , J.E. Myers , The Problem of Physico-Chemical Periodicity (Arnold & Co , London, 1926 )
7. A. Lotka , Z. Phys . Chem. 72 , 508 ( 1910 )
8. J. Hirniak , Z. Phys . Chem. 75 , 675 ( 1910 )
9. A.J. Lotka , Z. Phys . Chem. 80 , 159 ( 1912 )
10. D. Aubin , A.D. Dalmedico , Hist. Math. 29 , 273 ( 2002 )
11. K. Tomita , T. Kai , F. Hikami , Prog. Theor. Phys. 57 , 1159 ( 1977 )
12. C. Knudsen , J. Sturis , J.S. Thomsen , Phys. Rev. A 44 , 3503 ( 1991 )
13. K. Tomita , Phys. Rep. 86 , 113 ( 1982 )
14. B.L. Hao , Elementary Symbolic Dynamics and Chaos in Dissipative Dynamical Systems (World Scientific , Singapore, 1989 )
15. B.-L. Hao , W.-M. Zheng , Applied Symbolic Dynamics and Chaos (World Scientific, Singapore , 1989 )
16. B. Hao , W. Zheng , Applied Symbolic Dynamics and Chaos , 2nd revised edn. (Peking University Press, Beijing, 2014 )
17. J.A.C. Gallas, Adv. Atom. Mol. Opt. Phys . 65 , 127 ( 2016 )
18. L. Junges , J.A.C. Gallas , J. Opt. Soc. Am. B 33 , 373 ( 2016 )
19. J.G. Freire , R. Meucci , F.T. Arecchi , J.A.C. Gallas, Chaos 25 , 097607 ( 2015 )
20. L. Junges , J.A.C. Gallas, Opt. Comm. 285 , 4500 ( 2012 )
21. L. Junges , J.A.C. Gallas, New J. Phys. 17 , 053038 ( 2015 )
22. J.G. Freire , J.A.C. Gallas, Phys. Chem. Chem. Phys . 13 , 12191 ( 2011 )
23. J.G. Freire , J.A.C. Gallas, Phys. Lett. A 375 , 1097 ( 2011 )
24. M.R. Gallas , M.R. Gallas , J.A.C. Gallas, Eur. Phys. J. Special Topics 223 , 2131 ( 2014 )
25. A. Sack , J.G. Freire , E. Lindberg , T. Po¨schel , J.A.C. Gallas, Sci. Rep. 3 , 03350 ( 2013 )
26. E. Pugliese , R. Meucci , S. Euzzor , J.G. Freire , J.A.C. Gallas, Sci. Rep. 5 , 08447 ( 2015 )
27. J.G. Freire , R.J. Field , J.A.C. Gallas , J. Chem . Phys. 131 , 044105 ( 2009 )
28. M.A. Nascimento , J.A.C. Gallas , H. Varela , Phys. Chem. Chem. Phys . 13 , 441 ( 2011 )
29. J.G. Freire , C. Cabeza , A. Marti , T. Po¨schel , J.A.C. Gallas, Sci. Rep. 3 , 1958 ( 2013 )
30. J.G. Freire , C. Cabeza , A. Marti , T. Po¨schel , J.A.C. Gallas, Eur. Phys. J. Special Topics 223 , 2857 ( 2014 )
31. V. Kovanis , A. Gavrielides , J.A.C. Gallas Eur. Phys. J. D 58 , 181 ( 2010 )
32. M.J.B. Hauser , J.A.C. Gallas , J. Phys. Chem. Lett . 5 , 4187 ( 2014 )
33. M.R. Gallas , J.A.C. Gallas, Chaos 25 , 064603 ( 2015 )
34. R.J. Field , J.A.C. Gallas , D. Schuldberg , Comm. Nonlin. Sci. Num. Sim . 49 , 135 ( 2017 )