Early warning signals for critical transitions in a thermoacoustic system

Dynamical systems can undergo critical transitions where the system suddenly shifts from one stable state to another at a critical threshold called the tipping point. The decrease in recovery rate to equilibrium (critical slowing down) as the system approaches the tipping point can be used to identify the proximity to a critical transition. Several measures have been adopted to provide early indications of critical transitions that happen in a variety of complex systems. In this study, we use early warning indicators to predict subcritical Hopf bifurcation occurring in a thermoacoustic system by analyzing the observables from experiments and from a theoretical model. We find that the early warning measures perform as robust indicators in the presence and absence of external noise. Thus, we illustrate the applicability of these indicators in an engineering system depicting critical transitions.

established in the system could cause structural damage, reduce performance of the power generating system or even result in disastrous events such as operational failure of space rockets 31 . Thus, early detection of the transition to self-sustained oscillations is critical in a combustion system.
In general, physical systems display complex dynamics which makes it difficult to predict the transitions or even to describe the nature of the transitions observed in such systems. Therefore, controlled laboratory-scale experiments are necessary to verify the robustness of the early warning measures in physical systems. In this study, we employ ingenious laboratory-scale experiments along with a theoretical model to study the effectiveness of the early warning indicators. This is the first attempt to apply these early warning measures to predict subcritical Hopf bifurcation in a prototypical thermoacoustic system. The experimental set-up consists of a horizontal duct with a heat source located inside the duct. The experimental system could undergo a transition to self-sustained pressure oscillations for suitable choice of its control parameters. In addition to experiments, we use a theoretical model derived from the conservation equations of momentum and energy.
In this study, we identify early warning signals such as increase in variance and conditional heteroskedasticity close to a subcritical Hopf bifurcation. We find that these measures are robust even when the system under consideration is under the influence of external noise. Many physical systems undergo critical transition from a non-oscillatory state to an oscillatory state, which can be described via a subcritical Hopf bifurcation. Hence, our study on the viability of the early warning measures to predict transition in a physical system gains considerable significance.

Results
Experimental set-up. The experimental set-up consists of a 1 m long duct of 10 cm × 10 cm square cross-section (Fig. 1). The duct is open to atmosphere at one end. The other end is connected to a rectangular chamber of size much larger than the cross-section dimensions of the duct. This chamber is referred to as decoupler 32 . A blower in the suction mode is used to establish the air flow through the duct. The decoupler isolates the duct from the fluctuations upstream. Thus, the pressure at both the ends of the duct are equal to the ambient pressure. Therefore, the pressure fluctuations at the boundaries are negligible.
An electrically heated wire mesh located in the upstream half of the duct acts as a source of heat for this prototypical system, unlike a flame in practical combustion systems. The voltage drop across the mesh and the current flowing through the mesh are measured to estimate the electrical power supplied to the mesh. We used loudspeakers to provide external fluctuations in experiments when we studied the robustness of early warning measures in the presence of noise. Summary of the model. The theoretical model considered in this study captures the feedback between the fluctuating heat release rate and the sound waves (pressure and velocity fluctuations). This feedback could be modified, by varying certain control parameters, resulting in a state shift to large amplitude periodic oscillations. The formulation of the model is detailed in the Methods section. The following set of equations, derived from the conservation equations of momentum and energy, describes the temporal evolution of the system 33 .
The model described in equation (1) represents a nonlinear self-excited oscillator and it displays a transition from a non-oscillatory state to an oscillatory state via subcritical Hopf bifurcation for suitable change in control parameters. Although the experimental system under study is laminar, aperiodic fluctuations exist in the flow field. These fluctuations in the flow field arise from the blower. The flow noise generated from the blower does not completely decay in the decoupler. These small amplitude fluctuations convect downstream giving rise to small amplitude fluctuations in the flow field. Thus, the base state of the system can be viewed as a laminar flow field superimposed with the low amplitude fluctuations arising from the blower. We chose additive Gaussian white noise 34 to model these background fluctuations in the theoretical model as it captures the dynamics observed in experiments qualitatively. This qualitative similarity between the experimental results and the theoretical model, when the source of fluctuations is modeled as additive Gaussian white noise, is reported both in the case of a prototypical thermoacoustic system 34 and in the case of an industrial gas turbine combustor 35 . Once we include the noise, equation (1) becomes: where, ε represents strength of the additive Gaussian white noise φ(t). The noise term φ(t) in equation (2) models both background fluctuations in the flow field and the external fluctuations imposed by the loud speaker. The sources of both these fluctuations; the fluctuations due to flow noise from the blower and the fluctuations induced by the loudspeaker; are external to the system. Hence, we use additive noise to model both these fluctuations. The non-dimensional intensity of background fluctuations and external fluctuations is referred to as β. The procedure to obtain β is detailed in the Methods section. In order to simulate the experimental conditions, we maintained the strength of additive Gaussian white noise ε in the model such that the non-dimensional noise intensity β is the same as that in the experiments 34 .
Critical transition associated with a subcritical Hopf bifurcation. A bifurcation analysis was performed in order to characterize the nature of the transition observed in experiments. The control parameter in this study, the heater power (K), was increased in a quasi-steady manner and the time series of unsteady pressure pertaining to each value of the parameter was recorded. When K was increased above a threshold value A, the median of the peak acoustic pressure (P) registered an abrupt increase (see Fig. 2).
In Fig. 2, the abrupt increase in P corresponds to a transition of the system from a non-oscillatory to oscillatory state. It is evident from earlier studies that the oscillatory state corresponds to a state of limit cycle oscillations 32,34 . Once the system reaches the oscillatory state, the parameter K should be reduced to a threshold value further below A, to bring the system back to the non-oscillatory state. The threshold values at which the transitions occur are different for the forward path (when K is increased) and for the reverse path (when K is decreased). This difference in threshold values results in the presence of a bistable regime. The presence of the bistable regime indicates that the transition to the oscillatory state observed in the experimental system is subcritical in nature 36 . The subcritical nature of Hopf bifurcation prompts us to use early warning measures to predict the occurrence of this transition.
In experiments, we have observed a 5% uncertainty in the value of the Hopf point. The uncertainty arises because the experimental conditions such as the mean temperature and the speed of sound are not exactly repeatable. Thus, in a particular experiment, we can expect the onset of oscillations at any K within the range of approximately 380 ± 20 W.
Critical slowing down indicators at the onset of critical transition. We calculated variance, lag-1 autocorrelation and also conditional heteroskedasticity of the time series in search for forewarning signals before the system reaches the state of self-sustained oscillations. Unlike a laboratory experiment where the control parameter is varied in a quasi-static way, the control parameter must be varied continuously as in the case of real systems. This prompted us to perform experiments by varying the control parameter in a continuous manner and capture the time series depicting the continuous transition from non-oscillatory state to oscillatory state. We varied K by 2 W in every 20 s and recorded the corresponding time series of acoustic pressure. A plot that depicts the change in K with time is included as Supplementary Fig. S3. The metric based methods suggested by Dakos et al. 19 Figure 2. The bifurcation diagram obtained from experiments. The non-oscillatory state loses its stability at point A (Hopf point) and the system undergoes a transition from non-oscillatory state to large amplitude oscillatory state as the parameter K is increased. The system returns to the non-oscillatory state at point C as K is reduced. We can observe a bistable region ABCD, where the system can remain either in a stable oscillatory state or in a stable non-oscillatory state. The inset shows, the variation of any measure M with control parameter μ obtained from the normal form equation associated with subcritical Hopf bifurcation. The stable and unstable branches are colored blue and red respectively. are adopted to calculate the early warning measures. The time series prior to the transition alone is used for calculating the early warning measures as we require precursors for an impending transition.
In Fig. 3a, the amplitude of the oscillations starts to grow only at around t = 644 s which marks the occurrence of Hopf bifurcation (at K = 394 W). The effects of critical slowing down, which can be observed just before the occurrence of the subcritical Hopf bifurcation, will aid in detecting the impending transition. Therefore, we used the time series up to t = 640 s to calculate the early warning measures. We can observe a significant increase in the value of variance well before the transition, which clearly serves as an early warning signal (Fig. 3). However, the lag-1 autocorrelation decreases before the transition, rendering it as a less effective early warning measure (Fig. 3).
We extended our search for effective early warning measures for predicting a subcritical Hopf bifurcation also in the theoretical model given by equation (2). We varied the non-dimensional control parameter (k) as a function of time, given by k = 1.667 × 10 −3 t, and obtained the corresponding time series of unsteady pressure. Here also, we observed a transition to a state of limit cycle oscillations 38 . The lag-1 autocorrelation and variance were calculated using the same procedure that we followed for the time series obtained from experiments. In the model, we observe that the values of variance and lag-1 autocorrelation increase well before the transition (Fig. 4). Lag-1 autocorrelation shows a decreasing trend in experiments while we observe an increasing trend in the model. In experiments, the system has background aperiodic fluctuations, though we have not added any external noise in the system. The presence of fluctuations makes autocorrelation a less effective precursor as it is established that autocorrelation can increase or decrease in the presence of fluctuations 39 . In order to use autocorrelation as a precursor, we need to have multiple realizations of the transition, which may not be feasible when it is applied as an early warning measure in a real time system 36 . The observations from the model and experiments make it clear that variance is a more robust early warning measure as compared to lag-1 autocorrelation. Our observation of variance being a more robust early warning measure than lag-1 autocorrelation is completely in conformity with Ghanavati et al. 40 .
The robustness of the early warning indicators must be tested in the presence of fluctuations as practical combustion systems work in a turbulent environment. Therefore, we added external noise of non-dimensional intensity β = 0.2 in the system using loudspeakers 34 . Here also, we changed the control parameter by 2 W in every 20 s and recorded the time series of unsteady pressure. We observe that lag-1 autocorrelation shows an increasing trend close to the transition and variance shows a concurrent rise indicating the proximity to a critical transition (Fig. 5). Our experimental results prove that early warning indicators perform very well even in the presence of external fluctuations of intensity one order higher than that of the background fluctuations.
We further tested the robustness of early warning signals using the model (equation (2)). In the model, we applied Gaussian white noise such that the non-dimensional intensity is 0.2 to match the experimental conditions. The value of the Hopf point (k = 0.62) does not change since the noise in the model is additive in nature 41 . The parameter k reaches the value of Hopf point only at t = 372. We can notice in Fig. 6(a) that the amplitude starts to grow much before the Hopf point is reached. Therefore, the transition to limit cycle presented in Fig. 6(a) represents noise induced tipping in the system. In this case also, we observed a trend similar to that observed in experiments, where variance and lag-1 autocorrelation increase well before the transition (Fig. 6). Our results show that the early warning measures perform well in forewarning critical transitions irrespective of the presence of background and external fluctuations.
We then proceed to calculate conditional heteroskedasticity which was introduced as an early warning indicator by Seekell et al. 18 . Conditional heteroskedasticity is indicated by the persistence in the conditional variance of the error terms. We followed the procedure suggested by Seekel et al. 18 to estimate conditional heteroskedasticity (refer Methods section for details of the procedure). In this procedure, the time series is modeled as an autoregressive process and the residuals are obtained. We then estimate the conditional heteroskedasticity by examining and (c) lag-1 autocorrelation as the system approaches the critical transition. We observe increase in variance and autocorrelation well before the transition. We maintained the value of ε such that the non-dimensional noise intensity is 0.02 to match the experimental conditions. Note that the parameter k and the acoustic pressure that we calculate from the model are non-dimensional. (c) lag-1 autocorrelation as the system approaches the critical transition. We observe a rise in variance and autocorrelation before the transition. We have added external noise to the system such that the non-dimensional noise intensity is β = 0.2. the persistence in the conditional variance of the residuals. We adopt conditional heteroskedasticity as an additional measure along with variance and autocorrelation to early warn the transition.
We consider the time series data used in Figs 3 to 6 to calculate the conditional heteroskedasticity. We identified the significant number of tests where conditional heteroskedasticity is observed. We find that the cumulative number of significant tests (denoted by C) for conditional heteroskedasticity, applied to the time series, increases as the critical transition is approached (Fig. 7). The increase in C for all the cases considered in this study ascertains regime shift in the system 18 .

Discussion
We find that the critical slowing down indicators, e.g., variance and conditional heteroskedasticity, as early warning signals are able to predict catastrophic transitions in a physical system. These measures are robust even in the presence of high intensity noise. In the model, lag-1 autocorrelation shows an increasing trend in the presence of low and high intensity noise. However, in experiments, lag-1 autocorrelation shows a decreasing trend (in the  presence of background fluctuations) as the critical point is approached, though it shows an increasing trend in the presence of external noise imposed by a loud speaker. The lack of robustness of autocorrelation as a precursor is due to the fact that we are using a single sample path or a single realization 36 .
This is the first experimental study where we establish the effectiveness of critical slowing down indicators in early warning critical transition in a thermoacoustic system. Our findings are highly pertinent as many real systems exhibit critical transitions. By implementing the early warning measures in systems such as the one considered in the present study, we can prevent critical transitions in real systems. We wish to emphasize that the early warning indicators exploit the phenomenon of critical slowing down which happens near any bifurcation. Therefore, any transition which can be viewed as a bifurcation can be predicted using the early warning measures.
There are certain limitations to the approach we have adopted in this study. The early warning indicators considered in the current study may fail to detect the regime shifts that occur as a result of rapid change in the control parameter 42 . Apart from this, the early warning measures might not detect transitions induced by a finite amplitude perturbation in the system. The perturbation can occur at any arbitrary instant of time and can take the system to the basin of attraction of the stable limit cycle. Tipping in the system as a result of perturbations is not preceded by critical slowing down. Hence, the early warning indicators discussed in this study might not detect this transition. We provide an example from our experiments and mathematical model where early warning measures fail to detect a transition induced by a sudden finite amplitude perturbation (see Supplementary material).

Methods
Data acquisition. We use piezoelectric transducers (PCB piezotronics, PCB103B02) to record the variation in acoustic pressure. The sensitivity of the pressure transducer is 217.5 mV/kPa and the resolution is 0.2 Pa. The output from the pressure transducer is acquired using a data acquisition system (PCI 6221) through a signal conditioner. The transducer is mounted 325 mm from the end open to the atmosphere. We add external noise to the system through loud speakers (Ahuja AU 60) mounted on the duct 625 mm from the end open to the atmosphere 35 . The white noise signal is created in LabVIEW software and then input to the loud speakers through an amplifier. A DC power supply (TDK-Lambda, GEN 8-400) of range 0-8 V and 0-400 A is used to heat the wire mesh. The wire mesh is made from mild steel. The uncertainty in the heater power is 0.4 W. The accuracy of the flow meter is ± 0.1%. The decoupler has dimensions 1200 mm × 450 mm × 450 mm.
Theoretical model. The model is derived from the conservation equations of momentum and energy 33 .
In the model, we consider a horizontal duct with a concentrated heat source located inside the duct. An average value of temperature in the duct is adopted to model the sound propagation in the duct. We can write the conservation equations of momentum and energy (in one-dimension) as follows 43 : where, γ is the ratio of specific heats of air. The quantities with tilde are dimensional. The variables pressure ( p), density (ρ ), velocity ( u) and heat release rate ( Q) are decomposed into mean and fluctuating quantities. For example, density is expressed as ρ ρ ρ = + ′  . The fluctuating quantities are much small compared to the mean. For instance, the amplitude of the pressure fluctuations is less than 1% of the mean pressure in Rijke tube. It implies that the propagation of sound waves in the duct can be modeled as a linear phenomenon. Therefore, the linear terms involving the fluctuating variables alone are retained and higher order terms are neglected. In addition, the low Mach number approximation 44 is adopted to obtain the following equations: The fluctuating heat release rate per unit volume ( ′  Q ) is given by the following relation 45 : where, λ is the heat conductivity of air, C V is the specific heat at constant volume, L W is the effective wire length, ρ and T are the mean density and temperature of the flow respectively, T W is the temperature of the wire, d W is the wire diameter and S is the cross-sectional area of the duct. The unsteady heat release rate is a nonlinear function of the velocity fluctuations at the heater location. The heat release rate responds to the velocity fluctuations after a delay (τ ) which is incorporated in equation (7). The physical reason for this delay is the presence of the thermal and the hydrodynamic boundary layers around the wire. Due to the presence of the boundary layer, the heat where c 0 is the speed of sound, L is the length of the duct, u 0 is the mean speed of the flow and d is the diameter of the cylindrical heating wire.
Moreover, a similar delay is expected in our experimental set-up, as we have a mesh type heater (which can be viewed as an array of cylindrical wires). In a practical combustion system, this delay can be attributed to factors such as the convective time scale (time taken by the fuel air mixture to travel the distance between the fuel injection point and the flame location), evaporation time scale, mixing time scale and reaction time scale.
The location of the heat source inside the duct is denoted by  x f . Although the sound propagation is modeled as a linear phenomenon, the interaction between heat release rate and the sound waves is nonlinear. Therefore, equations (5) and (6) still represent a nonlinear system of equations as the expression for heat release rate is nonlinear 33 . Equations (5) and (6) are then expressed in terms of non-dimensional quantities by adopting the following scales for non-dimensionalization.
The quantities with tilde are dimensional. M is the mean flow Mach number and p is the steady state pressure. The following equations are obtained after non-dimensionalization.
The variables that govern the heat conduction from the wire are encapsulated into a single non-dimensional parameter denoted by k. In the analyses using the model, γ = k k M 2 / , is used as the control parameter analogous to the heater power in experiments. The non-dimensional acoustic pressure and acoustic velocity are represented by ′ p and ′ u . The acoustic pressure ( ′ p ) and the acoustic velocity ( ′ u ) are expressed in terms of the natural spatial modes of the duct as follows 33 : Here η j and η  j represent the time varying coefficients of j th natural mode. The expansions for ′ u and ′ p from equation (12) are substituted into equations (10) and (11) 32,47 . ω j = jπ represents the non-dimensional angular frequency of the j th mode. The non-dimensional frequency ω j is obtained by non-dimensionalizing the angular frequency of the j th natural mode of a duct open at both ends, ω π =  jc L 2 ( /2 ) j 0 . ω  j is non-dimensionalized with c 0 /L, where c 0 /L is the travel time of the sound waves in the duct. The value of the parameters used for computation are N (number of modes) = 10, dt (step size) = 0.01, ω j = jπ, c 1 = 0.1, c 2 = 0.06, x f = 0.25 and τ = 0.2.
The amplitude of fluctuations or the noise level in the system (both in experiments and in the theoretical model) is estimated by measuring the rms amplitude of the acoustic pressure when the system is in the non-oscillatory state. The noise amplitude (I) is non-dimensionalized by the amplitude of limit cycle attained at the Hopf point in the absence of external noise (P H ).
Here, β is the non-dimensional noise intensity.
Derivation of early warning signals. We performed the statistical analyses using the "Early Warning Signals Tool Box" (http://www.early-warning-signals.org/). We used the time series up to the impending critical transition to calculate the early warning indicators. For variance and autocorrelation, we calculated the temporal trend by estimating the nonparametric Kendall rank correlation (τ k ) 48 Kendall's τ k is a statistical tool used to measure the association between two measured quantities. Lag-1 autocorrelation is determined by the autocorrelation function (ACF) given below: where, P(t) is the value of the state variable at time t, and P m and σ 2 are the mean and variance of P(t) within the time frame considered. Variance is calculated as: where N is the number of observations.
We used Lagrange multiplier test to calculate conditional heteroskedasticity 49 . For this, we fit autoregressive model of selected order to the time series. The order is selected according to Akaike information criterion 50 which is a measure of the relative goodness of the fit. After fitting of model to the time series, we calculated squared residuals (ξ t 2 ) and regressed them one time step where α 0 and α i represent the regression coefficients. The residuals give the error variance relationship which shows the properties of conditional heteroskedasticity. We conducted chi square test to compare the values of R 2 (Lagrange multiplier test statistic) to a χ 2 distribution so as to identify the number of significant tests where conditional heteroskedasticity is observed.