Resonant Behavior in a Periodically Forced Non-Isothermal Oregonator

3 hours ago - For the current study we considered the three-variable Oregonator model with the temperature incorporated as a variable (not a parameter...
0 downloads 0 Views 2MB Size
Subscriber access provided by Nottingham Trent University

A: Kinetics, Dynamics, Photochemistry, and Excited States

Resonant Behavior in a Periodically Forced Non-Isothermal Oregonator David García-Selfa, Alberto P Munuzuri, Juan Pérez-Mercader, and David S.A. Simakov J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b05238 • Publication Date (Web): 23 Aug 2019 Downloaded from pubs.acs.org on August 29, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Resonant Behavior in a Periodically Forced Non-Isothermal Oregonator David García-Selfaa, Alberto P. Muñuzuria,b, Juan Pérez-Mercaderb,c and David S. A. Simakov*,d Group of Nonlinear Physics, Universidad de Santiago de Compostela, Campus Sur, 15782 Santiago de Compostela, Spain a

b Department

of Earth and Planetary Sciences, Origins of Life Initiative, Harvard University, Cambridge, Massachusetts 02138, United States c Santa

Fe Institute, Santa Fe, New Mexico 87501, United States

d Department

of Chemical Engineering, University of Waterloo, Waterloo, Ontario N2L 3G1,Canada

______________________________________________________________________________

ABSTRACT: Nonisothermal chemical oscillators are poorly studied systems because chemical oscillations are conventionally studied under isothermal conditions. Coupling chemical reactions with heat generation and removal in a nonisothermal oscillatory system can lead to a highly nontrivial nonlinear dynamic behavior. For the current study we considered the three-variable Oregonator model with the temperature incorporated as a variable (not a parameter) thus adding an energy balance to the set of equations. The effect of the temperature on reaction rates is included through the temperature-dependent reaction rate coefficients (Arrhenius law). To model continuous operation in a lab environment, the system was subjected to external forcing through the coolant temperature and infrared irradiation. By conducting numerical simulations and parametric studies, we have found that the system is capable of resonant behavior exhibiting induced oscillations. Our findings indicate that an external source of heat (e.g., via an infrared light emission diode) can be used to induce a Hopf bifurcation under resonant conditions in an experimental Belousov–Zhabotinsky reactor.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. INTRODUCTION Oscillatory phenomena are common in Nature and observed in a great variety of biological processes such as the Krebs cycle, the circadian clock, and the cell cycle.1 A classic example of a chemical oscillator is the Belousov-Zhabotinsky reaction,2-3 which has become a prime example for non-linear, out-of-equilibrium chemical dynamics. Natural oscillators are normally subject to both periodic and stochastic external perturbations in the values of environmental physical variables, such as light irradiation and temperature. Such fluctuations are expected to modulate, modify, or even determine and entrain the behavior of nonlinear systems.4-5 A notable example of such entrainment from Biology is temperature compensation in the circadian clock, which has been analyzed theoretically employing the Goodwin oscillator model.6-8 Another notable example of nonlinear dynamic behavior induced by external forcing is resonance. This phenomenon is recognized as the emergence of periodic behavior in a dynamical nonlinear system subject to a subthreshold periodic forcing, sometimes in the presence of stochastic modulation (stochastic resonance) or even in the absence of any periodic forcing, induced purely by noise (coherence resonance). Such phenomena have been observed in mathematical models4 and experimentally, using the photosensitive Belousov-Zhabotinsky reaction.5 The resonant behavior is typically recognized as the presence of a maximum in a plot of oscillation amplitude versus forcing frequency (for periodic forcing) or noise amplitude (for stochastic resonance). In this paper we explore the effect of temperature modulation using an idealized model of the Belousov-Zhabotinsky (BZ) reaction, the Oregonator,9 which we have suitably modified to include a heat balance with temperature considered as a dependent variable via dynamical reaction heat evolution. The BZ reaction was first proposed as an oversimplified description of the Krebs cycle. 2 ACS Paragon Plus Environment

Page 2 of 25

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Since its discovery2 and understanding,3 it has been considered as a model of chemical oscillatory behavior, even to the point where the observations within this context were used to predict and control systems as complicated as the cardiac muscle.10 Several types of external forcing of the BZ reaction via temperature modulation have been explored both experimentally and theoretically.11-14 It has also been proposed that the BZ reaction could have a temperature compensation mechanism,11 suggesting that it is possible to entrain the behavior of a chemical oscillator by periodic modulation of temperature. A common feature of these studies on the effects of temperature modulation11-14 is that the temperature is considered as an externally forced parameter, not a dependent variable. Therefore, the effect of heat evolution was not investigated explicitly in these studies. In reality, heat effects, such as reaction heat generation, can transform the dynamics of a chemical system very significantly, especially in concentrated systems. This phenomenon is well described for high temperature, heterogeneous systems that can exhibit self-induced oscillating temperature pulses, e.g., in the catalytic oxidation of carbon monoxide.15-17 To the best of our knowledge, the effect of dynamic heat evolution has not been studied yet explicitly (considering temperature as dependent variable, not a parameter) in low temperature, homogeneous oscillators such as the BZ reaction. In what follows, we will consider the Oregonator model of the BZ reaction operated in a flow reactor as a model system and we will include a transient heat balance with the reaction heat, heat exchange (via cooling), and external forcing (via infrared irradiation) terms, in order to reflect nonisothermal conditions. This full description of a chemical system allowed us to introduce periodic fluctuations in the reaction mixture temperature indirectly, and therefore mimicking realistic conditions. We will demonstrate that the external periodic modulation via infrared irradiation can

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 25

trigger a non-trivial oscillatory behavior, namely induced pseudo-deterministic oscillations and resonant behaviour.

2. METHODS Let us consider a lab-scale, well-mixed flow reactor under unsteady-state conditions (a transient continuously stirred tank reactor (CSTR)) a configuration that was recently reported in an experimental study.5 Mass balances for the BZ reaction carried out in this reactor are modeled by the three-variable, five-step Oregonator scheme18 that includes a flow term:19

 dX 2 2  dt  k1 H AY  k2 HXY  k3 HAX  2k4 X  k f X   dY  k1 H 2 AY  k2 HXY  k5 fBZ  k f X  dt   dZ  dt  2k3 HAX  k5 BZ  k f Z 

(1)

In eq 1, X, Y, Z are variable molar concentrations: 𝑋 = [𝐻𝐵𝑟𝑂2] (activator), 𝑌 = [𝐵𝑟 ― ] (inhibitor), 𝑍 = [𝑀𝑜𝑥] (catalyst). Other concentrations (𝐴 = [𝐵𝑟𝑂3― ], 𝐻 = [𝐻2𝑆𝑂4] and 𝐵 is the total concentration of organic species) are considered as constant parameters as in previous publications on the Oregonator model with flow terms.19 𝑓 is a stoichiometric factor, 𝑘𝑓 is the feed rate, and 𝑘𝑖 are reaction rate constants. Temperature dependence of the reaction rate constants

[ (

𝐸𝑖 1 𝑇

follows the Arrhenius law: 𝑘𝑖(𝑇) = 𝑘𝑖0𝑒𝑥𝑝 ― 𝑅

1

)]

― 𝑇0 , with 𝑘𝑖0 = 𝑘𝑖(𝑇0) being the rate constant

at a reference temperature 𝑇0. Heat evolution is modeled using the transient energy balance for an externally cooled CSTR,20 which is also subject to external forcing via infrared irradiation (which can be easily implemented in an experimental setup using a light-emitting diode (LED)5): 4 ACS Paragon Plus Environment

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

VR  wC pw

dT  VR  H Ri ri  Q& VR RIR dt i

(2)

In eq 2, the term ( ― 𝑉𝑅∑𝑖∆𝐻𝑅𝑖𝑟𝑖) represents heat generation as a result of chemical reactions, the heat exchange rate is given by Newton’s law of cooling, 𝑄 = 𝑈𝑆(𝑇𝐶 ―𝑇), and the term (𝑉𝑅𝑅𝐼𝑅) represents reaction mixture heating by infrared irradiation. Tc is the coolant (e.g., chilled water) temperature and U is the heat transfer coefficient. See Supporting Information for details. Combining eqs 1 and 2 and using Tyson’s scaling,21 this non-isothermal Oregonator model transforms to the following system of equations:

 dx  d  x(1  x)  y (q  x)   x   dy  fz  y (q  x)   y  d   dz  x  z   z  d  d    ( x, y, z ,  )  (1   )    d

(3)

2𝑘4

𝑘2

𝑘4𝑘5𝐵

In eq 3, the dimensionless variables 𝑥 = 𝑘3 𝐻 𝐴𝑋, 𝑦 = 𝑘3𝐴𝑌 and 𝑧 = (𝑘

3

𝑍 represent the

𝐻 𝐴)2

dimensionless concentrations of the activator, inhibitor and catalyst respectively, θ=T/Tc is the 𝑘𝑓

dimensionless temperature, and 𝜏 = 𝑘5𝐵𝑡 is the dimensionless time. 𝜅 = 𝑘5𝐵 is the dimensionless 𝑘5𝐵

2𝑘4𝑘5𝐵

flow rate and other dimensionless parameters are 𝜖 = 𝑘3𝐴𝐻, 𝛿 = 𝑘 𝑘 𝐴𝐻2, 𝑞 = 2 3

2𝑘1𝑘4 𝑘2𝑘3

and 𝛾 =

𝑉𝑅𝜌𝑤𝐶𝑝𝑤𝑘5𝐵 𝑈𝑆

. The term Ψ(x,y,z,θ) in the heat balance represents the reaction heat, the term (1-θ) is the heat exchange, and the term 𝜔 =

𝑉𝑅𝑅𝐼𝑅 𝑈𝑆𝑇𝑐

represents the external source of heat. See Supporting

Information for notation and values of parameters. Note that, since the reaction rate constants depend on temperature, all dimensionless 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

parameters in eq 3 (, , , q and ) are temperature-dependent in a non-linear way (via the Arrhenius law). We can readily see that there are two mechanisms to introduce perturbations in the temperature of the system. One mechanism is additive in character, via the irradiation intensity (). Another mechanism is multiplicative, via the coolant temperature (1-θ). In the present study, only additive perturbations will be considered.

3. RESULTS Numerical simulations were carried out using previously published activation energies and rate constants,22 which are listed in Supporting Information, alongside with other parameters. The reactor volume was set to 50 mL, representing a lab-scale reactor where excellent mixing and heat transfer can be achieved.5 The system of ordinary differential equations described by eq 3 was integrated using the 4th order Runge-Kutta method, using the following initial values: x = 0, y = 0.1, z = 0, θ = 1. The control parameters are the H2SO4 concentration in the solution (H0), the external coolant temperature (Tc), and the heat transfer coefficient (U). We first analyze the system behaviour in the absence of external forcing ( = 0 in eq 3) setting the heat transfer coefficient to U = 10 J m-2 K-1 s-1. For the H0 range used (0.07-0.16 M), Hopf bifurcation takes place for Tc = 283-313 K (10-40 C), Figure 1. For a given H0, the system switches to an oscillatory state when Tc passes a certain threshold value. For higher H0 values, Hopf bifurcation occurs at higher Tc. In Figure 1, the reactor temperature (T) was equal to the coolant temperature (Tc) after a short transient, which is due to the excellent heat transfer because of the small reactor volume and relatively high heat transfer coefficient.

6 ACS Paragon Plus Environment

Page 6 of 25

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 1. Phase diagram for the non-isothermal Oregonator with reaction heat and cooling terms (in the absence of external forcing by infrared irradiation). Black solid line marks the Hopf bifurcation boundary.

Figure 2 demonstrates the Hopf bifurcation for H0 = 0.10 M. As the coolant (and reactor) temperature is increased, the amplitude of the oscillations is decreased until the stationary (nonoscillatory) state is reached, with the supercritical Hopf bifurcation occurring at 296.3 K. The inset in the figure shows the variation of the oscillation period versus the coolant temperature (Tc). Note that since after the Hopf bifurcation oscillations disappear, the oscillation period is only plotted until the bifurcation point is reached. For temperatures less than T = 288 K, the system will stay in oscillatory state for this particular values of H0 (Figure 1), with both the oscillation amplitude and period further increasing, until water freezing point is reached, where oscillations are physically impossible. Recall that T = Tc, as explained above with regards to Figure 1.

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

Figure 2. Hopf bifurcation: The oscillation amplitude (xmax - xmin) is plotted vs. coolant temperature (Tc) for A0 = 0.24, B0 = 0.18 and a given value of H0 (H0 = 0.10), showing the Hopf bifurcation at 296.3 K. The inset shows the oscillation period vs. Tc in the oscillatory state.

3.1 Sinusoidal Forcing with Infrared Irradiation We now consider the same system (an externally cooled lab-scale CSTR) but subject to external forcing introduced, for example, via an infrared (IR) LED. A similar configuration was recently reported in an experimental study of the photosensitive BZ reaction.5 The intensity is set to vary periodically with time. This effect is introduced into the model via the external forcing parameter (𝜔 =

𝑉𝑅𝑅𝐼𝑅 𝑈𝑆𝑇𝑐

) in eq 3. The transient light intensity of the IR LED, 𝑅𝐼𝑅, is simulated according to the

following equation (RIR0 is the baseline, Af and ff are the forcing amplitude and frequency): 𝑅𝐼𝑅 = 𝑅𝐼𝑅0 + 𝐴𝑓 cos (2𝜋𝑓𝑓𝑡)

(4)

From the remainder of this section we set Tc = 283 K and H0 = 0.1 M; these values correspond to the oscillatory state in the absence of external forcing, Figure 1. However, when the external 8 ACS Paragon Plus Environment

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

forcing is applied the reactor temperature increases shifting the system towards the Hopf bifurcation boundary. Note that in Figure 1 the reactor temperature (T) was equal to the coolant temperature (Tc) after a short transient (due to the small reactor volume and excellent heat transfer). Therefore, Figure 1 also represents the phase diagram in the H0 vs. T plane. In our analysis below, the intensity of the IR LED (𝑅𝐼𝑅) is periodically varied forcing the temperature of the reactor to vary accordingly in such a way that the system periodically crosses the Hopf bifurcation. We considered U=10 J m-2 K-1 s-1 (as in Figure 1) and RIR0 = 0.028 kW/L, Af = 0.018 kW/L in eq 4. When RIR oscillates between 0.010 and 0.046 kW/L the Hopf bifurcation occurs at 296.3 K (Figure 2). The natural oscillation frequency in the vicinity of the Hopf bifurcation is f0 ≈ 0.00459 Hz. Figure 3 shows three cases corresponding to three different forcing frequencies: ff = 0.000045 Hz (case 1), ff = 0.00459 Hz (case 2), ff = 0.0135 Hz (case 3). The left column in Figure 3 shows the temporal evolution of the reactor temperature (as a consequence of the periodic variations of the IR intensity reaching the reactor) and the column on the right shows the evolution of the dimensionless activator concentration. Temporal evolution of the dimensionless concentrations of the inhibitor and catalyst is shown in Supporting Information, Figure S1 (for case 1). The reference cases with RIR0 = 0 kW/L (no forcing) and RIR0 = 0.018 kW/L (constant forcing) are shown in Figures S2 and S3 in Supplementary Information. Without forcing, the system is in the deterministic oscillatory state (Figure S2). Constant forcing initially results in the temperature increase (and decaying deterministic oscillations), until the balance between the heat supply (irradiation) and removal (cooling) is established, eventually leading to a stationary, nonoscillatory state at T = 297 K, which is above the Hopf bifurcation (Figure S3).

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Periodic (sinusoidal) forcing: Left column shows the evolution of the temperature in the reactor for three different values of the forcing frequencies (ff). The right column shows the evolution of the (dimensionless) activator concentration. Dashed lines in left column corresponding to the Hopf bifurcation (THB = 296.3 K). Parameters: H0 = 0.1 M, Tc = 283 K, U = 10 J m-2 K-1 s-1, RIR0 = 0.028 kW/L, Af = 0.018 kW/L. 10 ACS Paragon Plus Environment

Page 10 of 25

Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4. Resonant curves obtained for several values of the heat transfer coefficient (U), showing the reactor temperature (left column) and the dimensionless activator variation amplitude (right column) in the range of normalized forcing frequencies. The natural frequency is f0 = 4.5910-3 Hz. Parameters: H0 = 0.1 M, Tc = 283 K, RIR0 = 0.028 kW/L, Af = 0.018 kW/L.

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4 presents a systematic analysis in the range of the forcing frequencies normalized to the natural frequency near the Hopf bifurcation (ff /f0), for several values of the heat transfer coefficient (U). The left column in Figure 4 shows the maximum (circles) and minimum (stars) values of the reactor temperature as a function of the normalized forcing frequency. As expected, higher temperature is obtained as the magnitude of U decreases (lower heat transfer coefficient). The plots on the right column show the difference between the maximum and minimum values of the activator as a function of the forcing frequency. The upper branches in Figure 4b,d represent deterministic oscillations, which are obtained for low frequencies, when the Hopf bifurcation is crossed (such as in the insert in Figure 3b).

Figure 5. Resonant curves for several values of the heat transfer coefficient U. The inset in the figure shows the resonant forcing frequency versus the heat transfer coefficient. Parameters: H0 = 0.1 M, Tc = 283 K, RIR0 = 0.028 kW/L, Af = 0.018 kW/L.

12 ACS Paragon Plus Environment

Page 12 of 25

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Different rows in Figure 4 correspond to different values of the heat transfer coefficient (U). Decreasing the value of U results in a less pronounced resonant peak. Note that the temperature at the resonance peak becomes larger when U is decreased (less efficient heat transfer). This phenomenon is clearly seen in Figure 5 that compares resonance curves for different values of U. The inset in Figure 5 shows the resonant forcing frequencies versus the heat transfer coefficient.

4. DISCUSSION We now analyze the three different cases shown in Figure 3, with the forcing frequency being below, similar to, and above the natural frequency: cases 1, 2, and 3 respectively. Case 1 (upper row, forcing frequency below natural frequency) shows trivial behaviour: the system responds with oscillations when the Hopf bifurcation is crossed. High forcing frequency (last row, case 3) results in a stationary state after the dissipation of a transient. In this case, the system dynamics is not fast enough to respond to high frequency fluctuations and only senses the averaged IR irradiation intensity that leads to the reactor temperature above the Hopf bifurcation, thus a stationary state. The non-trivial behavior occurs for the intermediate frequency, which is close to the natural frequency near the Hopf bifurcation (second row, case 2). In this case, the system exhibits oscillations despite the fact that the reactor temperature is always above the Hopf bifurcation (Figure 3c) and, therefore, oscillations are expected to arrest. The maximum and minimum temperatures shown in Figure 4 strongly depend on the forcing frequency, decreasing/increasing for increasing frequencies, until overlap completely for higher frequencies. This behaviour reflects the finite heat capacity of the reaction medium (LHS of eq 2) that affects the rate of dynamic response of the system to external perturbations. When the reaction temperature is considered as a parameter, not a dependent variable, such behaviour cannot be

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

observed and the temperature can be apparently varied even at high frequencies, which is physically impossible. Therefore, considering the reaction temperature as a dependent variable allows for a more accurate representation of the system dynamics, especially at high forcing frequencies when the heat capacity effect becomes prominent. The plots on the right column show the difference between the maximum and minimum values of the activator as a function of the forcing frequency. For relatively low frequencies (ff /f0 < 0.1), the temperature oscillates between the two values and so does the activator variation amplitude. Deterministic oscillations, which are shown as upper branches in Figure 4b,d (red stars), only occur when the temperature is below the Hopf bifurcation (such as in Figure 3a,b). This behaviour is trivial. Above certain frequencies, oscillations are suppressed because the forcing frequency is much faster than the system dynamics. The forced variation amplitude of the activator (blue circles) is still not zero because the periodic changes in the reactor temperature force the system variables to adjust to the periodically changing conditions. The non-trivial behaviour is observed in the vicinity of the natural frequency (ff /f0  1), where the higher amplitude oscillations in the activator are observed (Figure 4b), although the reactor temperature does not cross the Hopf bifurcation point (Figure 4a). These oscillations have a relatively low amplitude, which is still significantly higher than the passive forced oscillations occurring at other frequencies, as it can be clearly seen from the sharp peak of the variation magnitude appearing at ff /f0  1. Note that for ff /f0  1 the reactor temperature does not oscillate and it is always above the Hopf bifurcation, in the non-oscillatory domain (Figure 4a). Therefore, the observed oscillations are not deterministic in their nature and this is clearly a resonant phenomenon that allows for the emergence of macroscopic oscillations. Note that the Oregonator model represents the Belousov-Zhabotinsky reaction, which is a nonlinear chemical oscillator. In 14 ACS Paragon Plus Environment

Page 14 of 25

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

addition, this oscillator is always damped in our model due to the cooling term (see eq 2). Therefore, the observed resonant behaviour is not expected to be identical to that observed in classical linear oscillators, where the disastrous resonance occurs at ff /f0 = 1 and in the absence of damping. Figure 5 shows resonant behaviour for different values of U. Each maximum of the variation amplitude curve (xmax – xmin) corresponds to a resonance. As in the classical forced oscillator with damping,23 resonances occur at lower normalized forcing frequency (ff /f0) and with lower amplitude for increasing damping (lower values of U), see the insert in Figure 5. Note that for lower values of U, cooling is less efficient and the reactor temperature increases, shifting the dynamic system away from the Hopf bifurcation towards the steady state region. This pattern coincides exactly with that of a damped forced oscillator for which the damping is increased.23

5. CONCLUSIONS We have considered the three-variable Oregonator model of the oscillatory Belousov-Zhabotinsky reaction in a flow reactor and, for the first time for this model system, incorporated the temperature as a dependent variable through the addition of the transient heat balance. The heat balance equation included the reaction heat generation and cooling term, as well as an additional term reflecting external forcing via infrared irradiation. The coupling between heat and material balances was included through the temperature-dependent reaction rate constants (via the Arrhenius law). The resulting four-variable non-isothermal Oregonator model not only accounts for changes in reaction temperature explicitly, in contrast to when forcing temperature is a parameter, but it also allows for investigating dynamic system response to fluctuations in external heat fluxes. 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Using experimentally obtained kinetic parameters for the Belousov-Zhabotinsky reaction from the literature and a set of representative parameters describing a lab-scale flow reactor, we were able to predict the occurrence of a Hopf bifurcation at realistic temperatures (10-40 C). We used additive perturbations to the system in order to control its dynamics. In that sense, periodic modulations of the infrared irradiation intensity were considered as a means of affecting the reactor temperature. The election of this control parameter is due to the simplicity of a potential experimental implementation, as such forcing can be simply done with an infrared light-emitting diode or laser. The results of numerical simulations have shown that, for some specific forcing frequencies, the system resonates and is able to exhibit sustained oscillations, even though the non-oscillatory state is the most probable configuration due to the parameters selected. Resonant behaviour was observed for forcing frequencies close to the natural oscillation frequency near the Hopf bifurcation. Higher order resonances were not observed due to the damping effect of the heat equation. It is, thus, possible to say that the periodicity of the environmental heat fluxes can induce periodic behavior in a chemical oscillator. This finding opens avenues for controlling chemical oscillators and opens a window for the extension to biochemical oscillators.

SUPPORTING INFORMATION Non-isothermal Oregonator model equations, values of the parameters used in simulations, typical values for the activator, inhibitor and catalyst, and simulation results for baseline cases (without forcing and with constant forcing).

AUTHOR INFORMATION Corresponding Author 16 ACS Paragon Plus Environment

Page 16 of 25

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

* Email: [email protected]. Phone: +1-519-888-4567.

ACKNOWLEDEMENTS The authors acknowledge funding support from the Natural Science and Engineering Research Council (NSERC) of Canada through the NSERC Discovery Grant program. APM and DGS are supported by the Spanish Ministerio de Economía y Competitividad and European Regional Development Fund under contract RTI2018-097063-B-I00 AEI/FEDER, UE, and by Xunta de Galicia under Research Grant No. 2018-PG082 and are part of the CRETUS Strategic Partnership (AGRUP2015/02) supported by Xunta de Galicia. All these programs are co-funded by FEDER (UE).

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

REFERENCES 1. Murray, J. D., Mathematical Biology. Springer: 1989. 2. Belousov, B. P., A periodic reaction and its mechanism. Ref. Radiats. Med. Medgiz, Moscow 1959, 145. 3. Zhabotinskii, A. M., Periodical oxidation of malonic acid in solution (A study of the Belousov reaction kinetics). Biofizika 1964, 9, 306-311. 4. Simakov, D. S.; Pérez-Mercader, J., Noise induced oscillations and coherence resonance ia a generic modelo f the nonisothermal chemical oscillator. Sci. Rep. 2013, 3, 2404. 5. Simakov, D. S. A.; Pérez-Mercader, J., Effect of noise correlation on noise-induced oscillation frequency in the photosensitive Belousov−Zhabotinsky reaction in a continuous stirred tank reactor. J. Phys. Chem. A 2013, 117 (51), 13999–14005. 6. Ruoff, P.; Vinsjevik, M.; Monnerjahn, C.; Rensing, L., The Goodwin oscillator: On the importance of degradation reactions in the circadian clock. J. Biol. Rhythms 1999, 14, 469–479. 7. Bodenstein, C.; Heiland, I.; Schuster, S., Temperature compensation and entrainment in circadian rhythms. Phys. Biol. 2012, 9, 036011. 8. Francois, P.; Despierre, N.; Siggia, E. D., Adaptive temperature compensation in circadian oscillations. PLoS Comput. Biol. 2012, 8, e1002585. 9. Field, R. J.; Noyes, R. M., Mechanisms of chemical oscillators: conceptual bases. Acc. Chem. Res. 1977, 10, 214–221. 10. Krinsky, V. I., Selforganization: autowaves and structures far from equilibrium. Springer: 1984. 11. Ruoff, P., Antagonistic balance in the Oregonator: about the possibility of temperaturecompensation in the Belousov-Zhabotinsky reaction. Physica D 1995, 84, 204–211. 12. Masia, M.; Marchettini, N.; Zambrano, V.; Rustici, M., Effect of temperature in a closed unstirred Belousov-Zhabotinsky system. Chem. Phys. Lett. 2001, 341, 285–291. 13. Bansagi, T. J.; Leda, M.; Toiya, M.; Zhabotinsky, A. M.; Epstein, I. R., High frequency oscillations in the Belousov-Zhabotinsky reaction. J. Phys. Chem. A 2009, 113, 5644–5648. 14. Novak, J.; Thompson, B. W.; Wilson, M. C. T.; Taylor, A. F.; Britton, M. M., Low frequency temperature forcing of chemical oscillations. Phys. Chem. Chem. Phys. 2011, 13, 12321–12327. 15. Flytzani-Stephanopoulos, M.; Schmidt, L. D.; Caretta, R., Steady state and transient oscillations in NH3 oxidation on Pt. J. Catal. 1980, 64, 346–355. 16. Sheintuch, M., Oscillatory states in the oxidation of carbon-monoxide on Platinum. AIChE J. 1981, 27, 20–25. 17. Liauw, M. A.; Somani, M.; Annamalai, J.; Luss, D., Oscillating temperature pulses during CO oxidation on a Pd/Al2O3 ring. AIChE J. 1997, 43, 1519–1528. 18. Field, R. J.; Noyes, R. M., Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction. J. Chem. Phys. 1974, 60, 1877. 19. Zhong, S.; Xin, H., Internal signal stochastic resonance in a modified flow Oregonator model driven by colored noise. J. Phys. Chem. A 2000, 104, 297-300. 20. Rawlings, J. B.; Ekerdt, J. G., Chemical Reactor Analysys and Design Fundamentals. Nob Hill Publishing: 2002. 21. Tyson, J. J., Oscillations, bistability, and echo waves in models of the Belousov-Zhabotinskii reaction. Ann. N. Y. Acad. Sci. 1979, 316, 279. 22. Pullela, S. R.; Cristancho, D.; He, P.; Luo, D.; Halla, K. R.; Cheng, Z., Temperature dependence of the Oregonator model for the Belousov-Zhabotinsky reaction. Phys. Chem. Chem. Phys. 2009, 11, 42364243. 23. French, A. P., Vibrations and Waves. The M.I.T. Introductory Physics Series. W. W. Norton & Company Inc.: New York 1971. 18 ACS Paragon Plus Environment

Page 18 of 25

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC GRAPHIC

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Phase diagram for the non-isothermal Oregonator with reaction heat and cooling terms (in the absence of external forcing by infrared irradiation). Black solid line marks the Hopf bifurcation boundary. 203x152mm (150 x 150 DPI)

ACS Paragon Plus Environment

Page 20 of 25

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Hopf bifurcation: The oscillation amplitude (xmax - xmin) is plotted vs. coolant temperature (Tc) for A0 = 0.24, B0 = 0.18 and a given value of H0 (H0 = 0.10), showing the Hopf bifurcation at 296.3 K. The inset shows the oscillation period vs. Tc in the oscillatory state. 242x182mm (150 x 150 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Periodic (sinusoidal) forcing: Left column shows the evolution of the temperature in the reactor for three different values of the forcing frequencies (ff). The right column shows the evolution of the (dimensionless) activator concentration. Dashed lines in left column corresponding to the Hopf bifurcation (THB = 296.3 K). Parameters: H0 = 0.1 M, Tc = 283 K, U = 10 J m-2 K-1 s-1, RIR0 = 0.028 kW/L, Af = 0.018 kW/L. 190x212mm (150 x 150 DPI)

ACS Paragon Plus Environment

Page 22 of 25

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Resonant curves obtained for several values of the heat transfer coefficient (U), showing the reactor temperature (left column) and the dimensionless activator variation amplitude (right column) in the range of normalized forcing frequencies. The natural frequency is f0 = 4.59×10-3 Hz. Parameters: H0 = 0.1 M, Tc = 283 K, RIR0 = 0.028 kW/L, Af = 0.018 kW/L. 190x219mm (150 x 150 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Resonant curves for several values of the heat transfer coefficient U. The inset in the figure shows the resonant forcing frequency versus the heat transfer coefficient. Parameters: H0 = 0.1 M, Tc = 283 K, RIR0 = 0.028 kW/L, Af = 0.018 kW/L. 199x181mm (150 x 150 DPI)

ACS Paragon Plus Environment

Page 24 of 25

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

217x162mm (150 x 150 DPI)

ACS Paragon Plus Environment