Periodic Motion in the Chaotic Phase of an Unstirred Ferroin-Catalyzed Belousov Zhabotinsky Reaction

The Belousov Zhabotinsky reaction, a self-organized oscillatory color-changing reaction, can show complex behavior when left unstirred in a cuvette environment. The most intriguing behavior is the transition from periodicity to chaos and back to periodicity as the system evolves in time. It was shown that this happens thanks due to the decoupling of reaction, diffusion and convection. We have recently discovered that, as the so-called chaotic transient takes place, periodic bulk motions in form of convective cells are created in the reaction solution. In this work we investigated this phenomenon experimentally by changing cuvette size and reaction volume, in order to allow different types of convection patterns to appear. So far, we have observed single and double convection cells in the system. There are indications that the convection patterns are connected to the duration of the chaotic phase. A simplified mathematical model confirms the form and dynamics of the observed convection cells and explains the connection between chemical chaos and hydrodynamical order.


OPTICAL SETUP
To follow the reaction dynamics a simple optical setup was used consisting of a ceiling LED panel (B.K.Licht 18 Watt LED Panel, 2400 lm, 4000K white light) as background light and a standard video camera (JVC Everio HM400) for data aquisition which was recorded on a PC (details see Fig. S1, (a)). Three cuvettes (two standard, one special cuvette (middle)) were put in the focal point of the video camera. The two standard cuvettes on each corner where inclinded slightly to adjust to the optical distortion of the video camera. The camera itself was placed on a solid support fixed with rubber bands. The dimensions of the used cuvettes are shown in Fig. S1, (b).

DATA TREATMENT
The video signal was recorded with the video software VLC media player with an aquisition rate of 0.5 to 1 images per second at high definition (HD) (1280 x 720 pixels). The images then were treated with the image processing package Fiji (imageJ2). The video camera JVC Everio HM400 uses a CMOS-Chip with Bayer RGB Color Filter Array. The blue filter in this array has a sensitivity peak at around 450 nm [9]. Ferriin, the oxidized blue form of the color indicator has a minimum at 450 nm [10]. This is the reason why we selected the blue channel of the video signal to analyze the dynamics of the reaction system. The information of a single channel is coded in 8-bit (256) different gray values from 0 (white) to 255 (black). High gray level values therefore correspond to low concentration of ferriin.

Extraction of the Time Series from the Video Data
To extract the time series from an video signal the video is transformed into a series of images. As mentioned above, the images were treated with the image processing package Fiji (imageJ2). This package is able to create a "Kymograph" (see Fig. S2 (b)), i.e. a space-time plot of the pixel values of a selected cross-section line (see Fig. S2 (a) green line) throughout an image series. In this way we have transformed the entire video into an image that contains the whole dynamics of the system at a certain line. In our case, it illustrates clearly the periodicity of the color change in the beginning and the gradually transition to chaos. Since we are interested in the dynamics of a given point (see Fig. S2 (a) red point) which corresponds to a line in the Kymograph (see Fig. S2 (b) red line) we translate the shades of gray back into gray level values by performing a "Profile Plot" and obtain Fig. S2 (c).
The time series extracted from the video data in the above described way is comparable to a time series obtained by a spectrophotometer at a wavelength of 450 nm. The advantage of an image series of the system over a spectrophotometric time series is that a time series can be in principle extracted at any point in the solution.

Visualization of Convection Cells
To visualize the convective motion in one image only, a time average over 60 seconds was calculated. The images were treated with the image processing package Fiji (imageJ2) where the option "image > stack > Z project" was used. Here as "Projection type" "Standard deviation" was used. In supplementary material S2 a video of the full reaction dynamics is shown (200x faster than real time).

Calculation of the Lyapunov Exponents
To ensure that the so called "chaotic" phase in our time series shows indeed typical properties of chaotic behavior we analyzed the corresponding time series thoroughly. A typical property of chaotic behavior is a broad band spectra in the frequency domain of the corresponding aperiodic oscillations as shown in the main article in Fig. 7 A-D. Another typical property of chaotic behavior is a positive maximal Lyapunov exponent. Rossi et al. have already shown that this is the case in the phase of aperiodic oscillations in ferroin-catalyzed unstirred BZ reactions [6]. To show that the chaotic phases in the time series obtained from the video data of our experiments shows also this property we calculated the maximal Lyapunov exponent. For this procedure the TISEAN Nonlinear Time Series Analysis package [3] was used. For calculating the Lyapunov exponent of the chaotic phase in the experiment shown in Fig. 1 of the main article (an experiment conducted in a fully filled standard cuvette), the first 1000 seconds of the chaotic phase (1400-2400s) were used. Fig S3 shows Figure S3. The local correlation entropy S(ε, n, m) as a function of the initial neighborhood size ε, the discrete time steps n and the embedding dimension M (for details see text) The fitted slope (blue dashed line in Fig. S3) for n in the range of (5,90) gives a value of 0.00256187 ± 0.0003971. This value is comparable to 0.00532 ± 0.00002 which is the value found by Rossi et al. for the chaotic phase in their experiment.
The above described procedure was repeated for the time points marked in the chaotic phase in the experiment that was conducted in the special cuvette (see Fig. 4 of the main article). All values for the maximum Lyapunov exponent were positive (λ A = 0.001178 ± 0.000276, λ B = 0.005099 ± 0.000249, λ C = 0.0009 ± 0.000271 and λ D = 0.001432 ± 0.000266). This confirms indeed the chaoticity of the chaotic phase.

THEORETICAL MODEL
In order to describe mathematically the above presented phenomena, we took into account the main processes which occur in the reaction cuvette. Firstly, we must consider the chemical reactions leading to localized inhomogeneities of the concentrations and the density. As a consequence, two other effects take place, namely the diffusion and convection of the liquid. In this model we neglect both the surface tension and the evaporation at the free upper surface of the liquid. The three processes taken into consideration are described by combining the classical differential equations for diffusion and reaction kinetics and the Navier-Stokes equation [4,2] which are given by: ∂ Ω ∂t Here c i are the concentrations of the representative reactants, . The corresponding chemical reactions rates f i (c) are given by Furthermore, u and Ω = ∇ × u represent, respectively, the fluid velocity and its vorticity, whereas P and ρ are the pressure and the density of the bulk, respectively. Because of the incompressibility of the bulk, the term ∇ u vanishes in equation S1b. The estimation of the Reynolds number, with respect to the whole cuvette, is about 5.52. Although this value is bigger than one, it is much smaller than typical values (of the order of 10 7 ). Consequently it is reasonable to neglect the nonlinear term in u in the Navier -Stokes equation (S1b) simplifies to Any way the solution obtained in formula (S24) leads to automatic disappearance of the nonlinear terms. Now, applying the Boussinesq approximation (adapted for the concentration instead of the temperature [5]), the concentrations slightly deviate from the homogeneous average and can be written as where c i denotes the average value of the concentrations and c i its deviations. Even if c i would be time dependent, we would have Also, because of the presence of these concentrations deviations, the bulk density may have small inhomogeneities and the pressure may deviate from its equilibrium value. Denoting with ρ the average density and with P 0 the pressure at the equilibrium, we have By keeping only the first order terms of the perturbing prime quantities, Eqs.(S1a) and (S1b) become Here, we are interested in performing a linear analysis around the steady state of equation (S7) instead of its full mathematical solution. Being at the steady state means that ∂ t c i = 0 and ∂ t Ω = 0. Also, thanks to the Boussinesq approximation, we have f i (c) = 0 (remember that c i are homogeneous). Consequently, we obtain where c * 3 remains an undetermined, free parameter. The chemical rates f i may be developed around c as follows Substituting in system (S7) and indicating with a star the quantities referring to the steady state, we get ν∆ where M i express the molar masses of the chemical species. To continue the calculations, we refer to the geometrical system depicted in Fig. S4.
x y z a b Figure S4. Vertical slab representing the cuvette used in the model. a and b being the dimensions of the liquid.

Two-dimensional Approximation
Since the width can be considered negligeable compared to the other two dimensions, we will consider that c 1 , c 2 , c 3 and u depend just on (y, z). Then, the dimensionless version of the system (S10) gets the form where u * τ = (u * 2 , u * 3 ) is the steady state velocity expressed in the dimensionless coordinates ξ = (ξ 2 , ξ 3 ), a ij are the dimensionless variante of ∂f i /∂c j | c and ω * 1 is the single non-vanishing component of the dimensionless vorticity Ω. The numbers Ra i represent the Rayleigh-numbers given by [5]: where g, L 0 , ρ, D i and ν are, respectively, the gravitational acceleration, the length gauge, the density, the diffusion coefficient and the kinematic viscosity. Four our system we have As the hydrodynamical model is reduced to two dimensions, it is sufficient to consider only the stream function Ψ, related to the velocities in the following manner The four unknown functions c 1 , c 2 , c 3 and Ψ are related by the following equations (hereafter we write c i for c i and u for u * ) It is interesting to note that ω 1 = −∆Ψ. The Laplace operator ∇ 2 Ψ is calculated only with respect to ξ 2 and ξ 3 . As it is known, the reaction rates f i (c) [see equation (S2)] are dimensionless quantities, so the involved constants are pure numbers = 10 −2 , q = 9.5 × 10 −5 , δ 1 = 1 , δ 2 = 0.558 , δ 3 = 0.959 .
The stoichiometric coefficient f ∈ (0.495, 2.4) for the periodic chemical regime [8]. The initial concentrations are homogeneous and have the following values The complete formal solution of the system (S15) will be addressed to a forthcoming publication. As mentioned above, in this article we intend to give qualitative hints explaining the onset of the convection cells and their relation to the periodic (or chaotic) chemical regime. The pattern of the convection cells is given by the velocity field (u 2 , u 3 ). Control parameters of the chaotic/periodic chemical regime may be obtained by the steady state values of c * 3 and the free stoichiometric coefficient f. Let us consider the velocities u τ = (∂ ξ 3 Ψ, −∂ ξ 2 Ψ) in equation (S15) as free unknown parametric functions and make the Ansatz c i = c 0i e κ· ξ (S18) in equation (S15a). This implies a secular equation for the paramater κ · u τ having the solution Though u τ is a given function on ξ, the r.h.s. of the foregoing expression shows us that k · u τ is independent of ξ. As a consequence, the term k · u τ can be considered as a number representing a constraint between the coordinates of the flow velocities. Depending on k · u τ , we can express the amplitudes of the deviations of the concentrations from the homogeneous average in the steady state when c 01 is still unknown. Substituting the previous result in equation (S15b), we obtain Remembering that ω * 1 = −∆Ψ, we obtain the equation for the stream function Ψ in the stationary regime which solution is The term Φ(κ, R) corresponds to The number R is determined by the following integral The former expression implies that R lies within the interval [min(a, b), The vertical walls of the cuvette interact similarly with the liquid, so we may inmpose the following boundary condition implying a "quantization" of κ 2 with n an odd integer (S28) The glass bottom of the cuvette and the free liquid surface are by no means physically equivalent. As stated above, we do not consider the superficial tension and the evaporation at the upper free liquid surface. Instead we impose that, at the two boundaries, the function Ψ behave in an "opposite manner" This implies also the "quantization" of the second component of κ with m an odd integer (S30) Of course this view of the boundary conditions is only one among many others given the multitude of other possible constrains. We chose this one because it corresponds to the observed dynamics in the cuvette used in this model. According to Saltzmann [7,1], we can consider a fraction of all possible modes of Ψ( ξ) The lowest mode is identified by n = m = 1. Consequently, we obtain The velocity field is easily obtained This solution corresponds to the observed direction of the flow (see Fig. S5 (b)). Considering the general solution (S24), we can easily verify that Substituting this expression in the secular equation (S19), we get a relation between the two control parameters c * 3 and f The foregoing expression enable us to qualitatively compare the hydrodynamic chaos with the chemical order. The smaller the numbers n and m are, the clearer is the structure of the convection pattern and the flow is more and more ordered. The variation interval of c * 3 is strongly affected by formula (S35). For fixed n and m, indeed, the range of c * 3 is determined by the possible variation of the stoichiometric coefficient f. Unfortunately, the fixed point expressed by equation (S8) does not admit, at least in the linear analysis, any oscillation around it, within the Marchettini model (S2). It is possible to show that FKN model is a limit case of Marchettini model if we take Scott [8] prooved that the fixpoint x * in FKN model has instability when f lies in the interval If we choose then we may find the instability range for c * 3 , including oscillations. Combining the formulae (S37) and (S38), we obtain the possible oscillation range c * 3 , which for q = 9.5 × 10 −5 ≈ 10 −4 gives We take this interval as a criterium to decide whether our system evolves to the chemical order or not. The threshold found by Marchettini et al. [5] is c 3 ≈ 0.93: under this level, the existence of a limit cycle has been noticed and above c 3 ≈ 0.95 the presence of quasi-periodicity and chaos. In our case, for n = m = 1 (hydrodynamic order) and for our dimensionless cuvette dimensions a = 1500 and b = 2000 c * 3 ∈ (6 × 10 −8 , 1.23 × 10 −7 ), thus outside of the chemical order. For n = 10 and m = 1, we are still outside: c * 3 ∈ (5.16 × 10 −6 , 1.17 × 10 −5 ) but for n = 100 and m = 1 we have c * 3 ∈ (5.6 × 10 −4 , 1.22 × 10 −3 ), thus already within the instability interval. These results prompt qualitatively the hypotheses that the hydrodynamic order is connected to chemical chaos, as seen from the experimental pictures (see Fig. S5).

Extension to Three Dimensions
The equation (S33) models only the flow dynamics in the special cuvette which is theoretically considered to be two-dimensional. The flow dynamics in the standard cuvette (see e.g. Fig. 5 of the main article) is three-dimensional, therefore a generalization of the model to three dimensions is needed. The corresponding generalization of the system (S11 a,b) has the form: is the steady state velocity expressed in the dimensionless coordinates ξ = (ξ 1 , ξ 2 , ξ 3 ). The same Ansatz (S18) for the concentrations leads to the same secular equation (S19) for κ u τ , having the same solution (S20) and consequently the same expressions (S21 a,b) for the concentration amplitudes c 0i . The involved quantity κ u τ has now the generalized meaning and we posses no stream function Ψ more. Instead we have the three conponents of the velocity obeying the following differential equations: The solutions of these equations are (we skip in the following the stars) given by where nowχ To regain the components of the velocity we apply the rotor operator on the definition of the vorticity ω = ∇ × u and use the fact that ∇ u = 0 for an incompressible fluid. So we find As ω is already known form (S43) solving the Poisson equation (S45) we obtain: Using the representation of the velocity field, and of the vorticity by means of their Fourier transforms one may find another relationship between u and ω: The remaining coordinate of the velocity u 3 may be found integrating the equation ω 2 = ∂u 1 ∂ξ 3 − ∂u 3 ∂ξ 1 : where g(ξ 2 , ξ 3 ) is an arbitrary function. These apperently redundant complications lead together with the fact that the fluid is incompressible to the results κ 1 = 0 and (S49) Through (S49 b) we refound the old result given by (S34). The condition that through the cuvette at least one flow circle takes place lead to the 'quantization' rules κ 2 = iπn a with n odd integer (S50 a) Coming back to the representation (S46) and selecting form e κ ξ just e κ ξ = e κ 1 ξ 1 sin(k 2 ξ 2 ) sin(k 3 ξ 3 ) we obtain (u 1 , u 2 , u 3 ) = 0,χ ( κ)φ( κ) k 2 k 3 sin(k 2 ξ 2 ) cos(k 3 ξ 3 ) 4π , −χ ( κ)φ( κ) k 2 2 cos(k 2 ξ 2 ) sin(k 3 ξ 3 ) 4π . (S52) Here k 2 = πn a , k 3 = πm b and the condition κ 1 = 0 was taken after performing the differential and integral computations in (S46). The appearing factorφ( κ) is given by: We stress the fact that the solution (S52) is only a sufficient one. It is implied by (S48) which itself is also only a sufficient, but not necessary condition. However it is able to describe the dynamichs in the standard cuvette at least in the case where it is partially filled.