1. Introduction

In recent years, there has been a noticeable rise in attention paid to investigating the intricate interplay between the complex dynamics and the chemical systems inherent in oscillating reactions. These studies have revealed various chaotic and non-equilibrium phenomena [1], including mixed mode oscillations [2], complex oscillations, bursting oscillations [3], bistability [4], intermittency, quasichaotic behaviour within the reactions [5] and coupled chemical oscillators [6] revealing intriguing phenomena in the chemical systems [7]. Autocatalytic reactions exert a profound in'uence on the stability and behaviour of systems, engendering complex dynamics, multiple stable states, and periodic oscillations. These phenomena collectively enrich the captivating realm of chemical kinetics [8] and the exploration of nonlinear reactions [9]. The Brusselator is a theoretical chemical model introduced by I. Prigogine and his collaborators [10] to represent autocatalytic chemical reactions with spatial oscillations. It resembles the well-known Belousov-Zhabotinsky (BZ) reaction [11], both exhibiting non-equilibrium behaviour with pattern formation and chemical oscillations. Over the decades, the Brusselator model has been extensively studied, shaping ideas on oscillations and pattern formation far from thermodynamic equilibrium. It’s a well-established system in the realm of non-equilibrium instabilities, with numerous research papers dedicated to its exploration, and the model has been subject to extensive investigation, encompassing limit cycles, Turing-Hopf bifurcation [12] and coupled systems [13]. Tomita and his collaborators [14] introduced an external forcing term to the Brusselator model, enabling the study of its response to external perturbations. This addition unveils intriguing and complex dynamics, including the emergence of multistability [15], hysteresis and vibrational resonance [16], controlling chaos [17], non-quantum chirality [18], control of a quasiperiodic route to chaos [19], etc.

Events in dynamical systems that suddenly occur and exhibit an unusual dynamical phenomenon are referred to as rare events or extreme events. Extreme events (EEs) refer to the sudden and random increase in the magnitude of one or more of the state variables of the dynamical system. They encompass a wide range of natural and human-made phenomena, including tsunamis, cyclones, rogue waves [20], earthquakes, droughts [21], chemical explosions, 'oods [22], stock market 'uctuations [23], pandemics [24] and certain transmissible diseases [25]. Despite the low probability of these rare events occurring, their resulting losses can be substantial. The prediction and analysis of extreme events occurring in real-world problems have not been fully mastered yet. Recent studies have delved into understanding the occurrence of such behaviour in dynamical systems [26, 27, 28, 29]. Extreme events arising in nonlinear dynamical systems mimic those observed in many physical systems, including ↓ber optics, nonlinear optics, photonics [30, 31], ↓nancial systems [32, 33], electronic oscillators [34, 35], mechanical oscillators [36], chemical oscillators [37] and neural models [38]. While studying nonlinear dynamical systems, the trajectory of a dynamical system typically follows a bounded attractor. Occasionally, they deviate signi↓cantly, causing a large magnitude of amplitude, spikes, or bursts due to instability in their state space. These rare occurrences resemble changes in the range of time series produced by the systems, which are often characterised by the large deviations from the nominal behaviour of a system [39]. To con↓rm these occurrences, speci↓c thresholds are determined by statistical analysis. The study of extreme events in dynamical systems is a growing area of research, as scientists and engineers are interested in understanding the mechanisms that lead to these events and developing methods for predicting them [40], predicting extreme events in dynamical systems is a challenging task, even in deterministic systems. There have been some attempts to ↓nd early warnings of extreme events [41, 42], but success in this area has been limited so far.

These rare events pose challenges due to limited data, making models more di cult. Developing available mathematical models can help us to understand the mechanisms behind the extreme events and improve our ability to predict and manage them. Conversely, the immediate hardware realisation of the nonlinear circuits capable of producing a variety of chaotic oscillations presents a substantial challenge for the future integration of chaos-based information systems. Notably, the exploration of techniques for generating diverse intricate chaotic oscillations, such as EEs, using uncomplicated electronic devices has attracted considerable theoretical and practical attention. The inquiry into creating the EEs with distinct electrical characteristics is captivating, given their unique nature as electrical signals [34, 35]. In continuation of the above, this study places its exclusive emphasis on the emergence of extreme events or occasional large-amplitude oscillations within the Brusselator chemical model when subjected to external periodic perturbations and investigates the dynamical instability causing these events. In the realm of chemical oscillators, only a limited number of experiments have provided evidence of their chaotic nature. To further investigate in the present manuscript the extreme events emerging from the Brusselator chemical model, we have designed a simple analogue electronic circuit and studied their dynamics experimentally in the laboratory. To the best of our knowledge, no prior research has explored this phenomenon within the context of this model or chemical oscillators.

This paper is organised as follows: In Section 2, we introduce the mathematical model of the proposed system and analyse its linear stability. In Section 3, we present the numerical study of the system by varying the system control parameters and also demonstrate the presence of extreme events in this model by representing unusually large events in time series using a qualifying threshold and occasional events in phase space. Further, we have plotted the probability distribution and return map to illustrate the rarity of these events as well as display the phase slips that occur during the extreme events. The global behaviour of the model is studied through the two-parameter bifurcation diagrams, which are given in Section 4. The experimental evidence of extreme events in this model is presented in Section 5. Finally, in sections 6, we draw our conclusions.

2. Mathematical Model of Brusselator

The Brusselator chemical model [10] is a trivial autonomous theoretical oscillator that mimics the autocatalytic reaction. The proposed autonomous Brusselator model can be written as,

where a and b are constants of system parameters and x and y represent the system state variables of the model (Eq. (1)). Then we have introduced an external periodic forcing term [14] to the Brusselator autonomous model of Eq. (1). The resultant system is a second-order nonautonomous di↑erential equation as follows:

The system of Eq. (2) is subjected to an external driving force with an amplitude f and frequency ω. Thus, the present model is composed of two simple ordinary di↑erential equations with external forcing, and the resulting oscillations are purely temporal (time-dependent). For understanding the dynamics of our proposed system, we integrated Eq. (2) using the fourth-order Runge-Kutta method with integration time steps dt = 0.01, initial condition x(0) = 0.01, and y(0) = 0.02.

2.1. Linear Stability Analysis

The stability criterion and the nature of bifurcational analysis on the Brusselator model were reported [10, 43]. The equilibrium point for the system Eq. 2 is obtained by setting the derivatives of x and y with respect to time equal to zero and without any external force. The eigenvalues of the system are λ_1,2 = (-(x² + b + 1 - 2xy) ). Thus, the system of Eq. 2 has only one equilibrium point, which is (x₀, y₀) = (a, b/a). The eigenvalues at the equilibrium point are λ_{1,2 equ} = (-(a² - b + 1) ). The stability of the system depends on the relationship between b and a² + 1. If b > (a² + 1), then , this indicates that the equilibrium point is unstable. If b < (a² + 1), then , this indicates that the equilibrium point is stable. If b = (a² + 1), then . This indicates that the equilibrium point is either a centre or an elliptic point, with oscillatory behaviour.

The stability plot for the (a - b) plane for λ_{1,2 equ} is shown in Fig. 1. The di↑erent colour zones are represented by di↑erent stability states, say unstable node (USN), unstable focus (USF), stable focus (SF), and stable node (SN). The critical line that represents centre stability is mentioned in an orange line in between USF and SF. From this stability analysis, the system stability changes with respect to parameters a and b since there exists a Hopf bifurcation for b = a² + 1.

3. Numerical Results

To begin with, numerically integrate and analyze the extreme events (EEs) in a system of Eq. (2) by examining its dynamics through a bifurcation diagram depicted in Fig. 2(a) and its corresponding maximum Lyapunov exponent in Fig. 2(b) by varying the system control parameter b within the range of b ∈(1.0, 1.2) and chosen the other system parameters are constant as a =0.2, f =0.06, and ω =0.7 with the initial conditions as x₀,y₀ = (0.01, 0.02). The bifurcation diagram was computed to detect the large peak values of state variable x in each value of the parameter b scanning. To distinguish the extreme events (EEs) from the normal events, a threshold height (H_s) of large amplitude oscillation is determined statistically [44], This statistical approach enables the identi↓cation and analysis of the occurrence of EEs in the system. The threshold value H_s is obtained by considering the dynamical aspects of the system observable, which involves measuring the signi↓cant large deviations away from the mean value of the state variable of the system. In other words, the EEs are identi↓ed as having large amplitudes substantially higher than the average value. The H_s is calculated as follows: H_s =< x_max > +nσ_xmax, where < x_max > represents the mean of the peak values of state variable x, σ_xmax is the standard deviation of the x_max, and n is an integer speci↓c to the system, typically ranging from 4 to 8 for extreme events, whose value determines and distinguishes extreme events from bounded chaos. To calculate H_s for a long run with iterations of 10⁸ time units, the system dynamics must be allowed to evolve past transient states. In this Fig. 2(a), we plotted the threshold height H_s for n =4 and n =6. By setting a threshold height, the extreme events can be distinguished from the normal events based on their large amplitudes, with EEs being those peaks that exceed the H_s value. In this one-parameter bifurcation diagram shown in Fig. 2(a), the system undergoes a transition from a regular attractor to a chaotic attractor via the usual period doubling (PD) route. The system’s dynamical summary, derived from the bifurcation analysis and accompanied by its corresponding Lyapunov exponents, is presented in Table 1 and shows the ranges within which the system’s dynamics exist.

The numerically obtained typical time series (x(t)) is shown in Fig. 3(a), and the corresponding phase portrait in the (x - y) plane of the bounded chaotic attractor is depicted in Fig. 3(b) for the value of b =1.1092. In Fig. 3(a) and 3(b), we observe that the chaotic behaviour remains con↓ned within a speci↓c region, characterised by relatively small amplitudes. Throughout the entire iteration, the system’s trajectory shows no evidence of surpassing the threshold value, where we set the basic criterion for extreme events (EEs) as n =4. The threshold value of the bounded chaos for the state variable x(t), calculated from H_s given earlier, is H_s(n = 4) =0.4112 for b =1.1092. Further, in the one-parameter bifurcation diagram of Fig. 2(a), when the control parameter b is varied, we ↓nd the occurrence of extreme events embedded in the bounded chaotic attractor. Interestingly, when the parameter b reaches the critical value of 1.1093, there is a sudden chaotic expansion observed in the x_max variable. This expansion is substantial enough to surpass the qualifying threshold (H_s) value. This phenomenon has been graphically depicted in Fig. 2(a) (with indicate arrow as EEs) for two distinct scenarios when n =4 (green colour) and n =6 (red colour). The typical extreme events can be clearly seen in the time series plot of the x(t) variable shown in Fig. 3(c) and the corresponding phase portrait in the (x - y) plane (Fig. 3(d)) for the value of b =1.1125. Here, the peaks in the variable x(t) with amplitudes four times and six times as large as the normal amplitudes are observed. In Figs. 3(c) and (d), we observe a signi↓cant expansion in the region of bounded chaos. This expansion in the time series exceeds calculated threshold values: H_s =1.7763 for n =4 and H_s =2.5012 for n =6, as indicated by the horizontal dashed red colour line in Fig. 3(c), and this sudden expansion in the system is characterised as an extreme events. As the control parameter b is varied further, the system dynamics bifurcate as non-extreme events with large amplitude chaotic behaviour are indicated as LAC in Fig. 2(a).

3.1. Statistical Characterization

We applied the statistical properties to characterise and con↓rm the extreme events (EEs) present in the Brusselator system of Eq. (2). In accordance with the ↓ndings in references [28, 32, 45], we have plotted the probability distribution function for the state variable x(t), we have taken the t-span length for a long run with iterations of 2 × 10⁷ time units, allowed the system to evolve past transient states for constant parameters a =0.2, ω =0.7, f =0.06 and various values of the parameter b and con↓rmed that the shape of the distribution does not change with respect to the t-span length. Fig. 4(a) shows the PDF for nominal chaos is bounded within a low range of x_max values below the H_s mark (vertical dashed black colour line), as expected for b =1.1092, signifying a non-extreme event. Fig. 4(d) reveals a continuously heavy tail distribution, and it surpasses the threshold represented in vertical dotted black colour line. This indicates the occurrence of extreme events, which con↓rms the low probability of the occurrence of large-amplitude events beyond the H_s mark (vertical black colour line) for b =1.1125, which is calculated for n =4 and n =6 using the respective temporal data.

Further, we have calculated the return map as an additional tool for di↑erentiating the bounded chaotic (BC) oscillations from the extreme events (EEs) oscillations. In Fig. 4(b),(e), we plotted the Poincar return map, which is obtained by plotting the ↓rst peak values x_maxn against the next peak values x_maxn+1 for the system (2). We observe a distinct behaviour in the bounded chaos for b =1.1092, where the return map remains locked and does not cross the H_s (plotted as a dashed red line). This con↓nement to the bounded region signi↓es the presence of non-extreme events, as shown in Fig. 4(b). Interestingly, as we change the parameter towards b =1.1125, we notice a signi↓cant increase in amplitude of x_max, leading to large events that cause the return map points to exceed the H_s (exceeding points are plotted as red hollow circles and H_s as a dashed red line). This crucial observation validates the existence of extreme events in system Eq. 2, as shown in Fig. 4(e). In addition, we plotted a projection of 3D phase plots in Fig. 4(c) and (f) to illustrate how the system transitions from bounded chaos (Fig. 4(c)) to extreme events (Fig. 4(f)) by including the external periodic force sin(ωt) as the third axis, and Fig. 4(f) represents the events that cross the threshold H_s indicated in red colour.

To better understand the system’s behaviour, we applied a Hilbert transform to the variable x(t) and an external periodic signal sin(ωt) as the third axis, as Fig. 4(f) represents. This allowed us to determine their instantaneous phase and subsequently ↓nd the di↑erence between these instantaneous phases, thus calculating the phase di↑erence δϕ. During the bounded chaos, we observed that the system remained phase-locked to the external periodic force. However, during the occurrence of large events, particularly extreme events, phase slips were found to occur because the system temporarily lost its phase-locking to the external periodic force, leading to a sudden and transient shift in its phase, which is shown in Fig. 5.

4. Two parameter bifurcation

The overall system dynamics across the parameter space (a - b) were depicted in Fig. 6, where and . Our focus is to identify and distinguish the regions of extreme events (EEs) from non-extreme events (NEEs). To begin with, we found that the EEs are exhibited when the system demonstrates chaotic behaviour. Hence, to distinguish the EEs region from the other region as a function of the parameters a and b, we estimated the threshold height (H_s). For a long run with iterations of 10⁸ time units, if the system shows nominal chaos with the maximum values of the x-variable (x_max) exceeding a threshold height for n =4, then the dynamics in that speci↓c region are marked as EEs (represented by the red colour points in Fig. 6(a). Conversely, if the system does not surpass H_s, then these parameter values are marked as part of the non-extreme events region (represented by grey colour points in Fig. 6(a). Throughout this analysis, we veri↓ed the emergence of these extreme events within the speci↓ed parameter ranges of a and b using time series data and one-parameter bifurcation diagrams. To con↓rm this with a two-parameter bifurcation for x_max depicted in Fig. 6(b), the dark region indicates a low height amplitude of the state variable x_max, suggesting the presence of periodic oscillations and bounded chaos in those regions. On the other hand, the colour gradient from violet to yellow represents high amplitudes of x_max. These high amplitudes indicate the occurrence of expansion in x_max. From Figs. 6(a) and 2(a), we can deduce that extreme events are found to be present only at the initial stage of the sudden expansion of x_max. As time progresses, the amplitudes tend to fall more frequently, leading them to no longer qualify as extreme events. Figures 6(a) and 6(b), it is evident that in the region where parameter and , there is a noticeable presence of high amplitudes (red-yellow colour) of x_max and a high probability of extreme events occurring in this range.

5. Experimental observations of Extreme Events

Recently, experimental circuit implementation has o↑ered an alternative avenue to investigate or verify the feasibility of physical implementations of theoretical dynamical models and to apply them in practical scenarios. To analyse the circuit dynamics and compare them with the numerical studies of a normalised model, it is derived by performing circuit variable substitutions and parameter transformations [35]. Therefore, in this section, we have designed a simple analog circuit and implemented it in the laboratory to validate the obtained numerical results of the two-dimensional nonautonomous Brusselator model given in the Eq. (2) described earlier.

According to the above numerical analysis, Fig. 7 shows the schematic of experimental circuit realisation in the top panel, and the complete analog circuit assembled using readily available discrete components was utilised to build the proposed chaotic circuit of the Brusselator model of Eq. (2), which is shown on a breadboard in the bottom panel. The resulting circuit is straightforward, cost-e↑ective, and can also serve laboratory experiments and educational purposes for exploring the innovative e↑ects in the dynamics of complex classical oscillators, as mentioned earlier. This circuit consists of linear resistors, capacitors, and μA741C operational ampli↓ers. The nonlinear functions in Eq. (2) are generated by using AD633JN analog multipliers. The operational ampli↓ers and multipliers operate with supply voltages of 12 V and saturated voltages of 9.5V. The arbitrary waveform generator (Agilent 33500B) is taken as the external periodic forcing voltage source f (t) = Fsin(Ωt).

By applying Kirchho↑’s circuit laws to the designed circuit (top panel) of Fig. 7, we get the following circuital equations:

Here, v_C1 and v_C2 are the voltages developed across the capacitors C₁ and C₂ respectively, represented by the circuit state variables.

Similar to the numerical simulation studies in section 3, we have observed the occurrence of EEs and bounded chaos followed by the period doubling route in experimentally also and the values of various circuit components are pre-determined in the circuit depicted in Fig. 7. For this breadboard experiment, speci↓c circuit component values can be chosen using an appropriate time scale [34]. The circuit’s time-constant-related circuit elements have been optimized to R = 100 kΩ and C = 2.2 nF. Further, we have ↓xed the values of other circuit elements as follows: The capacitance values of the capacitors are ↓xed as C₁ = C₂ =2.2 nF. The resistances of resistors are set as follows: R₁,R₃,R₆,R₇, R₈ = 10 kΩ, R₂ =6.8 kΩ, R₅ =3.3 kΩ, and R₄ = 100 kΩ. The amplitude F =4.5 V and frequency Ω =6.2 kHz of the external periodic voltage source. The gains of the analog multipliers M₁ and M₂ are (1∕10) V. It’s worth noting that all the circuit components have tolerances of 5 %. Also, we varied the resistance values of R₄ and R₆ in the breadboard circuit during the experimental observations and to match the selected control parameters a and b, as speci↓ed in Eq. (2), based on the results of numerical simulation studies.

Speci↓cally, in the laboratory experiments, a precision potentiometer is utilized for the adjustable feedback resistor R₄ in the ↓rst integrator (IC1) of the circuit (Fig. 7 (top panel), which is considered the control parameter. The tunable potentiometer has a resistance of 100 kΩ to adjust the external bias (E) in the experimental observations. The choice of these values is justi↓ed by our wish to use the same sets of system parameters for both numerical and experimental studies. The resistance of the variable resistor R₄ is gradually tuned; for di↑erent values of R₄, the circuit generates di↑erent attractors experimentally.

In order to begin our experimental study, in Fig. 7, when the resistance R₄ is set to 22.5 kΩ, the oscillator generates periodic limit cycle oscillations. When the control parameter R₄ is varied in the range of (23 kΩ ≤ R₄ ≤ 55.5 kΩ), we observed the period-doubling sequences; the range (55.5 kΩ ≤ R₄ ≤ 78.0 kΩ), chaotic regime, and the range (78.0 kΩ ≤ R₄ ≤ 82.6 kΩ), extreme events oscillations. The experimentally observed results are summarized in Figs. 8. Snapshots of the experimental time series and phase portraits are captured using the analog oscilloscope for di↑erent values of the circuit control parameter, the variable resistance R₄. The experimental circuit outputs can directly be displayed on an oscilloscope by feeding the output voltages v_C1 and v_C2 to connect them to the X and Y channels of the dual-channel analog oscilloscope (ScientiFic SM 410) with 1 V/div in the X direction and 1 V/div in the Y direction. The experimentally obtained time series of voltage v_C1(t) (Fig. 8a(i)) and corresponding the typical phase portrait of bounded chaotic attractor in the (v_C1 - v_C2) plane (Fig. 8a(ii)) for R₄ =56.5 kΩ. Continuing to vary the resistance R₄ in the range of R₄ (78.0 kΩ, 83.0 kΩ), we observe the occurrence of extreme events embedded within the bounded chaotic attractor, as shown in Fig. 8(b). The temporal time evolution plot for the voltage v_C1(t) and corresponding the typical phase portrait of extreme events in the (v_C1 - v_C2) plane for R₄ =80.5 kΩ as shown in Fig. 8b(i) and b(ii). Obviously, the experimentally captured results of chaotic behaviour and extreme events oscillations agree with those shown in Figs. 3 by the numerical simulations.

6. Conclusion

In this article, we have presented the dynamics of a Brusselator chemical system driven by an external periodic force, studied both numerically and experimentally. We have investigated and con↓rmed many dynamical phenomena, including the well-known period doubling, chaotic oscillations, and the most intriguing extreme events behaviour in the proposed chemical model. The one-parameter bifurcation diagram, maximum Lyapunov exponent, phase portraits, and time series segments characterise the system’s dynamics, while the probability distribution function, the instantaneous phase calculation, and the Poincar return map characterise the extreme events behaviour. Additionally, a two-parameter bifurcation diagram is used to distinguish extreme events within a two-parameter space. Further, the real-time hardware analog electronic experiments were carried out in the laboratory, and the ↓ndings con↓rmed the numerically obtained results. To the best of our knowledge, this is the ↓rst analog electronic experimental study on the Brusselator chemical model with the observation of extreme events. These results are giving more insight into how to construct real-time experiments with chemical reaction setups and avoid external stimuli. The system creates intermittency peaks [46] and multistability [47]. However, in our work, we have reported that external periodic forces cause very large oscillations and lead to extreme events for appropriate system parameter regimes. Adding the external periodic force to the Brusselator mathematical model is an initial assumption, and it helps to understand the real chemical experiments for the Brusselator autocatalytic reaction. The external periodic force has been applied to the state variable x. The autocatalytic nature of the system is evident in the state variables x and y, which evolve chaotically with extreme event oscillations. However, due to the addition of the external periodic force, the dimensions are increased, resulting in chaotic oscillations with rich dynamical behaviors. These intriguing results, such as extreme events and reverse period doubling, will be replicated in real Brusselator experiments to uncover the hidden dynamics in the near future.

Acknowledgment

K.T. acknowledges the Chennai Institute of Technology, Chennai, India, for its experimental and computational assistance. SS acknowledges the Basic Research Program of the National Research University, Higher School of Economics, Moscow.

Con'ict of Interest

The authors declare that they have no con'ict of interest.

DATA AVAILABILITY

The data that supports the ↓ndings of this study is available from the corresponding author upon request.

References