Effect of rate of change of parameter on early warning signals for critical transitions
EEffect of rate of change of parameter on early warning signals for criticaltransitions
Induja Pavithran and R. I. Sujith a) Department of Physics, IIT Madras, Chennai-600036, India Department of Aerospace Engineering, IIT Madras, Chennai-600036, India (Dated: 29 January 2021)
Many dynamical systems exhibit abrupt transitions or tipping as the control parameter is varied. In scenarios wherethe parameter is varied continuously, the rate of change of control parameter greatly affects the performance of earlywarning signals (EWS) for such critical transitions. We study the impact of variation of the control parameter witha finite rate on the performance of EWS for critical transitions in a thermoacoustic system (a horizontal Rijke tube)exhibiting subcritical Hopf bifurcation. There is a growing interest in developing early warning signals for tippingin real systems. Firstly, we explore the efficacy of early warning signals based on critical slowing down and fractalcharacteristics. From this study, lag-1 autocorrelation (AC) and Hurst exponent ( H ) are found to be good measuresto predict the transition well-before the tipping point. The warning time, obtained using AC and H , reduces withan increase in the rate of change of the control parameter following an inverse power law relation. Hence, for veryfast rates, the warning time may be too short to perform any control action. Furthermore, we report the observationof a hyperexponential scaling relation between the AC and the variance of fluctuations during such dynamic Hopfbifurcation. We construct a theoretical model for noisy Hopf bifurcation wherein the control parameter is continuouslyvaried at different rates to study the effect of rate of change of parameter on EWS. Similar results, including thehyperexponential scaling, are observed in the model as well. Critical transitions, which can result in sudden devastat-ing changes to the state of the system, are ubiquitous innatural, economic, and social systems. Continuous varia-tion of control parameter in time can delay the transitiondue to memory effects. Therefore, the rate of change ofcontrol parameter as well as the value of the control pa-rameter determine the tipping point. In such cases, weobserve a rate-dependent tipping-delay from the criticalpoint predicted by bifurcation analysis. Furthermore, theinterplay between the inherent noise in the system and therate of change of parameter blurs the stability margins.Therefore, it is essential to develop effective early warn-ings signals (EWS) to predict such transitions in practi-cal systems. In recent years, many studies have exploredvarious EWS to predict the onset of critical transitions.We conduct experiments in a thermoacoustic system ex-hibiting Hopf bifurcation. In a thermoacoustic system, thepositive feedback between the unsteady heat release ratefluctuations and the acoustic field in the confinement canresult in a transition to high amplitude limit cycle oscil-lations. These self-sustained large amplitude oscillationsare dangerous to the system. For different rates of changeof the control parameter, we compare the efficacy of var-ious early warning signals based on critical slowing downand fractal characteristics of the signals acquired from athermoacoustic system. We also investigate the variationof warning time with the rate of change of the parameter. a) Electronic mail: [email protected]
I. INTRODUCTION
The dynamics of many natural and human-made systemsare controlled by various parameters which evolve in time. Atiny perturbation in a system parameter can qualitatively alterthe state of the system when it crosses a critical threshold.This phenomenon is generally known as tipping or criticaltransition, wherein a small change of a parameter can cause asudden transition in the state of the system . The term tippinghas been used to explain various phenomena such as phasetransitions and bifurcations, which are often associated withdangerous and catastrophic transitions. Tipping is observed inreal systems such as climate systems , ecological systems , fi-nancial markets and biological systems . In earth’s climatesystem, a gradual change in local climate can affect ecosys-tems and can sometimes trigger a drastic switch to a contrast-ing state . Contagion in financial markets , spontaneousasthma attacks , and epileptic seizures are other instancesof tipping.There are various mechanisms through which tipping oc-curs. Recently, Ashwin et al . classified the underlyingmechanisms as bifurcation-induced (B), noise-induced (N) orrate-induced (R) tipping. B-tipping occurs when a system pa-rameter is varied slowly through the bifurcation point, com-monly known as the slow passage through the bifurcation. Ex-amples of B-tipping include saddle-node, transcritical, pitch-fork and Hopf bifurcations. In N-tipping, the system can jumpto another stable state due to the presence of noise of sufficientamplitude. In other words, noise can drive the system betweenthe coexisting attractors in systems exhibiting multistability.Sometimes, the rate of change of parameter plays a morepivotal role than the actual value of the parameter. Ashwin etal. showed that when a rate-sensitive parameter is varied asa function of time, at a slow rate, the system dynamics followsthe quasi-static attractor. For faster rates of change of the pa- a r X i v : . [ n li n . PS ] J a n rameter, above a critical rate, they observed that the systemcan be driven outside the basin of attraction of the quasi-staticattractor, and can evolve towards a new stable state resultingin rate induced tipping (R-tipping). On the other hand, byvarying the bifurcation parameter in a bistable system, onecan achieve preconditioned rate induced tipping, as demon-strated by Tony et al. . They reported that the system couldbe driven towards the basin of attraction of limit cycle beforethe actual loss of stability of fixed point, for fast enough rateswith a finite amplitude initial perturbation. Here, the tippingdepends on the rate of change of control parameter and ini-tial conditions. In these cases, the rate at which the parameteris varied determines the tipping point, not the absolute valueof parameter. However, practical systems may exhibit tippingphenomena as a result of a combination of bifurcation, rateand noise.In the current study, we explore the effects of rates ofchange of bifurcation parameter on ‘bifurcation induced tip-ping’. When we vary the bifurcation parameter continuouslyat a finite rate, tipping will be delayed considerably from theparameter value predicted by the bifurcation analysis . Dueto the continuous variation of the control parameter, the sys-tem stays in the vicinity of the unstable attractor for some timeeven after the stability is lost. This phenomenon is referred toas ‘rate-delayed tipping’ . The delay observed in the transi-tion is found to be dependent on the rate of change of param-eter and the initial conditions . Later, Majumdar et al. reported that this delay in tipping is independent of the rateof change of control parameter and determined solely by theinitial value of the parameter. Recently, Bonciolini et al. showed experimentally that the delay in bifurcation increaseswith rate of change of parameter. After such contradictoryobservations in literature, Unni et al. highlighted the role ofinterplay between the inherent noise in the system and rateof change of parameter in deciding the tipping point. Eventhough the delay increases with the rate of change of param-eter, noise brings a high variability in the tipping point. De-termining the stability margin is difficult for practical systemswhere stability boundaries are smeared due to this interplaybetween noise and rate. Devising efficient EWS for abrupttransitions in real systems is of critical importance. For exam-ple, predicting earthquakes, climate changes, and catastrophicevents in engineering systems are desirable from social andeconomic viewpoints. However, predicting such tipping be-fore the occurrence of the event is challenging because thesystem may not show any indication before the tipping pointis reached, especially when there is combined effects of rateand noise.Abundant studies discussing quasi-static bifurcations or B-tipping are available in literature. Despite the inherent char-acteristics of the systems, the dynamics close to the bifurca-tion point are found to be the same across different systems .The transitions through various types of bifurcations are re-lated and generic early warning signals exist for catastrophicbifurcations . The phenomenon known as critical slowingdown near the bifurcation point gives information about theimpending tipping for many types of bifurcations. As the sys-tem approaches the critical point, it recovers slowly from the external perturbations. This slow recovery leads to an increasein the memory of the system. The two commonly used earlywarning indicators that work based on critical slowing downare the lag-1 autocorrelation and the variance of the fluctua-tions of a system variable. These measures have been provento predict B-tipping, wherever the tipping is accompanied bya change of stability of the system .Recently, Wilkat et al. showed that there is no evidence ofcritical slowing down prior to human epileptic seizures. Theyconjecture that the tipping mechanisms for the human epilep-tic brain may be a combination of B-tipping, R-tipping andN-tipping and there may be no easily-identifiable EWS forsuch cases. Most tipping events occurring in nature involvesystem parameters changing at a constant, varying or undeter-mined rate along with considerable intensity of noise in thesystem . Then, tipping can be influenced either by noiseor by the effect of rate of change of the parameter. Thus,prediction becomes hard with conventional EWS. There aremany other issues when dealing with rate dependent phenom-ena. For example, consider the case where the entire transitionhappens at a very fast rate such that there is not enough dataavailable. Computing precursors in such conditions will bechallenging. Even if we get warning, will there be enoughtime to perform a control action? The outstanding question is:‘can we predict transitions in the real systems, considering thecombined effects of inherent noise and the rate of change ofthe parameter?’Many studies focused on continuous variation of pa-rameters employs numerical analysis of standard bifur-cation models and limited experimental studies areavailable . In the present study, we choose to workwith a prototypical thermoacoustic system (known as a hori-zontal Rijke tube) exhibiting Hopf bifurcation, because (i) weobserve a catastrophic transition similar to that observed inmany practical situations, (ii) we can obtain time series datafor a long duration with high sampling frequency, (iii) we canvary the control parameter at different rates and, (iv) we canrepeat the experiments at same conditions to verify the obser-vations.The Rijke tube undergoes a subcritical Hopf bifurcationfrom a non-oscillatory to an oscillatory state as we varythe control parameter . The oscillatory state with self-sustained large amplitude periodic oscillations is termed ther-moacoustic instability (TAI) in literature . In a confined ther-moacoustic system, the positive feedback between the acous-tic field and the unsteady heat release rate manifests as highamplitude acoustic pressure oscillations during TAI. The on-set of TAI can cause failure of rocket engines and damagegas turbine engines . Often, control parameters in practi-cal systems are changed continuously at a finite non-zero rate.Developing EWS for transition to TAI will help to evade suchdisastrous events. In this paper, we study the effects of rateof change of control parameter on the performance of variousEWS, by investigating the variation of warning time providedby EWS with the rate of change of parameter. Further, wepresent a mathematical model that captures qualitatively, thekey features observed in the experiments.The rest of the paper is organized as follows. Section II de- FIG. 1. Schematic of the horizontal Rijke tube used for the experi-ments. It comprises a 1 m long duct, an electrically heated wire meshand a rectangular chamber called decoupler. The wire mesh is shownseparately outside the duct. A piezoelectric transducer is mounted onthe duct to acquire the acoustic pressure fluctuations inside the duct. scribes the experimental setup. Subsequently, the results anddiscussions are detailed in Sec. III and the main conclusionsare given in Sec. IV. We provide the details of the method-ologies for calculating different measures in Appendix A. Weshow the robustness of EWS with the selection of thresholdfor warning in Appendix B and the analysis to check for falsewarnings in Appendix C.
II. EXPERIMENTS
We perform experiments on a laminar thermoacoustic sys-tem known as the horizontal Rijke tube (Fig. 1). The setupconsists of a horizontal duct with a square cross-section whichhouses an electrically heated wire mesh. Air enters the ductthrough a rectangular chamber known as the decoupler, whichisolates the duct from the upstream fluctuations. The decou-pler ensures that the pressure fluctuations at the end connectedto it are negligible. The duct is open to the atmosphere atthe end away from decoupler. Thus, the pressure at both theends becomes equal to the atmospheric pressure. This designhelps to maintain an acoustically open boundary condition(i.e., acoustic pressure fluctuations, p (cid:48) = 0) at both the ends. ADC power supply (TDK-Lambda, GEN 8-400, 0-8 V, 0-400A) is used to provide electric power to heat the wire mesh.The mass flow rate of air through the duct is controlled usingan electronic mass flow controller (Alicat Scientific, MCR se-ries) with an uncertainty of ± (0.8% of reading + 0.2% of fullscale). We measure the acoustic pressure fluctuations insidethe duct using a piezoelectric sensor (PCB103B02, sensitiv-ity: 217.5 mV/kPa, resolution: 0.2 Pa and uncertainty: 0.15Pa) at a sampling rate of 10 kHz. A more detailed descriptionof the setup can be found in Gopalakrishnan & Sujith .In the present study, we control the voltage ( V ) appliedacross the wire mesh and the current through the meshchanges accordingly. Therefore, we estimate the power ( P )generated in the wire mesh by measuring both the voltage and current. The uncertainty in the measurement of voltage is (0.1 V rated + V measured )% and the uncertainty in the measure-ment of current is (0.3 I rated + I measured )%, where V rated = 8 V and I rated = 400 A. All the other parameters such asthe location of the heater mesh (27.5 cm from the decoupler)and the mass flow rate of air (100 SLPM) are kept constantsuitably to obtain subcritical Hopf bifurcation.First, we perform experiments where V is varied in a quasi-steady manner by allowing the system to evolve for a finitetime duration at a constant P . We let the system reach itsasymptotic state and then measure the acoustic pressure fluc-tuations at different values of P . We select the maximumheater power to be less than 1152 W for all the experiments,as the wire mesh may melt and get damaged due to overheat-ing at higher powers. Subsequently, we increase V continu-ously at different rates and record the pressure signals duringthe ramp. Here, ramp refers to the continuous increase of theheater power in time. In this paper, we report a linear varia-tion of V ; i.e , the rate ( r = dV / dt ) of change of V is constant( V = V + r t ). As we have programmed V to vary linearly, I changes accordingly and then the corresponding variation of P is quadratic. III. RESULTSA. Rate-dependent tipping-delay in a thermoacoustic system
We first conduct quasi-static experiments to identify therange of parameter values required to capture the transitionto high amplitude limit cycle oscillations. We calculate theroot mean square ( rms ) value of the acoustic pressure fluctu-ations ( p (cid:48) rms ) acquired at different values of heater power ( P ).Figure 2a represents the bifurcation diagram showing the vari-ation of p (cid:48) rms as a function of P . We observe that p (cid:48) rms ≈ P corresponding to a quiescent state with amplitudelevels comparable to the noise floor ( ∼ p (cid:48) rms as an abrupt rise and is attributed to subcriticalHopf bifurcation. The parameter value at which this transitionoccurs, marked as µ in Fig. 2a, is known as the Hopf point.As mentioned earlier, we perform experiments with a lin-ear variation of V . Therefore, P changes continuously startingfrom a heater power which is far lower than µ and increasesthrough the Hopf point to a high value. Here, the variation of P is from 0 to 1152 W in a nonlinear fashion. Thus, through-out this study, we mention the rate of change of voltage withtime ( r = dV / dt ) which is kept constant for each run of theexperiment. In Fig. 2b-c, we plot the typical variation of V and P and the corresponding pressure signal ( p (cid:48) ) depicting thetransition from a quiescent state to high amplitude limit cycleoscillations.Figure 3 provides the evidence for rate-dependent tipping-delay as reported in literature . We plot the time evo-lution of p (cid:48) as a function of time-varying P for three different r (30 mV/s, 60 mV/s, and 120 mV/s) in Fig. 3. It is quite chal-lenging to define a tipping point for the onset of oscillations FIG. 2. (a) The variation of rms of the acoustic pressure fluctuations ( p (cid:48) rms ) for the quasi-static experiments. p (cid:48) rms shows an abrupt jump atthe heater power ∼
600 W, denoted as µ . (b) A typical variation of voltage ( V ) at a constant rate of increase ( dV / dt = 60 mV/s) and thecorresponding variation of heater power ( P ). (c) The pressure signal ( p (cid:48) ) acquired continuously as the control parameter is varied in time.FIG. 3. Time evolution of p (cid:48) as a function of time-varying P forthree different rates of change of V . The rate of change of rms offluctuations is also shown in the inset to identify the onset of TAI.Here, the maximum rate of change of p (cid:48) rms is considered as the onsetof oscillations. The delay in the transition, δ , is found to increasewith an increase in r . Here, p (cid:48) rms is calculated for a moving windowof 1 s with an overlap of 0.9 s. for dynamic bifurcations, unlike the quasi-static bifurcations.The difficultly in defining a tipping point is because the os-cillations take a finite time to grow, and the parameter wouldhave changed to another value by that time. Hence, we adoptthe following method to select the tipping point. We calculatethe rate of change of p (cid:48) rms . The sudden increase in the growthrate of oscillations reflects as a maximum in the rate of changeof p (cid:48) rms (shown in the inset of Fig. 3). The delay ( δ ) in the tip-ping is marked from the Hopf point to the maximum of rateof change of p (cid:48) rms (Fig. 3). Henceforth, we use this method todefine the onset of TAI in this study.According to bifurcation theory, the system loses its stabil-ity at the Hopf point, µ . Due to memory effects, the systemcontinues to be in the vicinity of the unstable attractor for afinite time. This phenomenon of rate-dependent tipping-delay occurs during slow passage through Hopf bifurcation as de-scribed by Baer et al. . On comparing the delay ( δ ) in theonset of TAI for different r , we observe that δ increases withan increase in r (see Fig. 3), congruent with the observationsreported in other systems . In the case of r = 120 mV/s(the fastest shown in Fig. 3), we observe a delayed onset of ∼
470 W from µ . Furthermore, we vary the control parameterat very fast rates, but limiting the values of P to 1152 W so asto not damage the heated wire mesh. In such cases, we do notobserve tipping within the duration during which the poweris varied; instead, it occurs when we let the system evolve atthe final value of P , allowing more time. Hence, we confirmthat the delay in the tipping increases with rate, by performingexperiments at various rates. This trend need not be the samefor highly turbulent systems where the inherent fluctuationscan perturb the unstable fixed point, causing it to tip towardsthe stable limit cycle. B. EWS for critical transitions
We compute several EWS for the pressure signal acquiredcontinuously during the ramp for the experiments performedat different rates of change of control parameter. The follow-ing measures are considered in this study: p (cid:48) rms , lag-1 auto-correlation (AC), variance (VAR), skewness (SKEW), kurto-sis ( K ), and Hurst exponent ( H ).Autocorrelation calculates the correlation of a signal withits delayed copy as a function of the delay. Lag-1 autocorrela-tion is the correlation between values in the signal that are onetime step apart. In this study, we refer to lag-1 autocorrelationas AC. Variance (VAR) is the expectation of the squared de-viation from the mean. For critical transitions, both AC andVAR increase based on the phenomenon of critical slowingdown. Systems approaching a transition to a new stable state,where the current stable state becomes unstable due to changein the control parameter, show slow response to external per-turbations. The phenomenon of slow recovery rate to the ex- FIG. 4. Variation of p (cid:48) rms , lag-1 autocorrelation (AC), variance (VAR), skewness (SKEW), kurtosis ( K ) and Hurst exponent ( H ) during thetransition to TAI in Rijke tube. Each column corresponds to the results for a particular rate ( r ) of change of V ((a)-(f): r = r = r = µ is marked on the x -axis. p (cid:48) rms , and VAR starts to increase almost at the onset of TAI, whereas, K and SKEW detect the tipping slightly before p (cid:48) rms increases. Incontrast, AC and H give early warning well-before the transition to TAI for all the cases shown here. ternal perturbations close to a critical transition is known ascritical slowing down (CSD). This slowing down leads to anincrease in autocorrelation and variance of fluctuations .When a system is driven gradually closer to a critical tran-sition, the increase in the autocorrelation occurs much beforethe actual transition. Similarly, the impact of perturbationsdo not decay fast, and their accumulating effect increases thevariance of fluctuations. CSD results in increased short termmemory of the system, and hence autocorrelation at low lagswould increase. CSD is observed in realistic models of spa-tially complex systems as well as in simple models andhas been used as EWS for critical transitions .Skewness (SKEW) and kurtosis ( K ) are not directly con-nected to critical slowing down. SKEW is a measure of thesymmetry of the probability distribution of the data about itsmean; i.e ., whether the distribution is biased towards one sideover the other. A distribution is said to be negatively skewed(or left skewed) when a majority of the data falls to the right of the mean. On the other hand, a distribution is positivelyskewed (or right skewed) when more data is concentrated tothe left of the mean. SKEW is zero for a normal distribution.Generally, asymmetry of fluctuations may increase (SKEWincreases) on approaching a catastrophic bifurcation, as thepotential landscape near the transition would be less steep onone side of the equilibrium than the other . Kurtosis ( K )gives information about whether the distribution has heavytails, or is more centred. The kurtosis of a normal distributionis 3, and a higher K indicates more outliers in the data. Veryclose to the tipping point, a comparatively longer distributionwith fatter tails results in a high value of kurtosis, K >
3. Adetailed procedure for the computation of these EWS is givenin Appendix A-3.Apart from these conventional measures, we investigate thefractal characteristics of the data close to tipping . We usea measure known as the Hurst exponent ( H ), which is re-lated to the fractal dimension ( D ) of the time series as H = − D . Fractal objects exhibit self-similar features atvarious scales of magnification; therefore, measures such aslength, area, and volume are dependent on the scale of mea-surement. For a fractal time series, the scaling of the rms ofthe standard deviation of fluctuations with the length of thedata segment gives H . For non-fractal objects such as si-nusoidal signals, H ≈ H need tobe selected carefully. As the system has preferred time scales(corresponding to the natural frequency) in the present study,selecting time scales less than one acoustic cycle may not cap-ture the periodicity present in the data and more than 4 acous-tic cycles may average out the fluctuations. Hence, we choose2 to 4 acoustic cycles for the calculation of Hurst exponent .We calculate H following the algorithm of Multifractal De-trended Fluctuation Analysis (MFDFA) , which is describedin detail in Appendix A-4. H is a measure of persistence or correlation of a signal. Ifan increase in the value is more likely to be followed by an-other increase in value, then the signal is called persistent. Apersistent or a positively correlated signal has H > .
5, ananti-persistent or negatively correlated signal (an increase invalue is mostly followed by a decrease in value, or vice versa)have H < . H = .
5. Frac-tal analysis has found a variety of applications in life sciences,engineering, econophysics, and geophysics . For in-stance, it has been used to distinguish healthy patients frompatients with heart failures . Similarly, the variations inthe H of geoelectric and seismic fluctuations provide indica-tors for earthquakes . In econophysics, Qian et al . used H as a measure of financial market predictability. In the presentstudy, we use H to predict critical transitions in a thermoa-coustic system.To study the effect of rate of change of parameter on theperformance of the aforementioned EWS, we plot the varia-tion of these measures as a function of P for three represen-tative cases of r : 5 mV/s, 30 mV/s and 80 mV/s in the first,second and third column, respectively (Fig. 4). We computethe measures for a moving window of size 1 s with an overlapof 0.9 s. The choice of this particular window size is to ensurethat at least 100 cycles of oscillations are covered in a window.The Hopf point ( µ ) is marked at 600 W with a red colouredarrow for reference, even though the tipping occurs after adelay. p (cid:48) rms indicates the growth in the amplitude of oscilla-tions during the transition. VAR detects the transition almostat the same P where the amplitude rises which is reflected asa steep increase in the p (cid:48) rms at the onset of TAI, while SKEWand K seem to perform slightly better than p (cid:48) rms and VAR.Initially, we observe small negative SKEW values indicatinga slightly left/negatively skewed distribution, and during thetransition, it shifts to a slightly right/positively skewed distri-bution. Nevertheless, there is no significant change in SKEWduring the transition, since − . < SKEW < . FIG. 5. Rate of change of p (cid:48) rms , lag-1 autocorrelation (AC), and Hurstexponent ( H ) with P is plotted for different rates of change of V .The maximum rate of change of AC and H occurs much before themaximum growth of p (cid:48) rms . As the rate of change of input voltagewith time ( dV / dt ) increases (from (a)-(d)), the delay (in terms of P )in growth of amplitude increases. considered as symmetric distribution. Thus, the distributionhas not changed in terms of skewness, and we cannot considerthe change in SKEW as a precursor for the tipping. Besides,kurtosis has a value ∼ K as a goodEWS for critical transitions in practical systems.Among these conventional measures, lag-1 autocorrelation(AC) appears to be the best EWS for this type of transition.An initial quiescent state results in near zero AC due to thelow amplitude uncorrelated noise. AC increases as periodic-ity increases and approaches 1 for limit cycle oscillations. Itis clear from these experiments that AC is a more robust earlywarning measure compared to p (cid:48) rms , VAR, SKEW and K . Ear-lier studies have reported that autocorrelation in the pres-ence of fluctuations is a less effective precursor compared tovariance. Actually, variance starts to build up gradually longbefore the transition, but it increases rapidly only at the tip-ping point. In fact, variance is the square of p (cid:48) rms . As we donot know the amplitude of oscillations in the final state, it isdifficult to rely on variance to determine when the transitionwill take place. In contrast, the value of AC is bounded be-tween 0 and 1, and we know how close we are to the tippingfrom the value of AC.The fractal based measure, H fluctuates between the val-ues 0.2 to 0.5 for the fixed point state. Even though the flowfield is laminar, there is noise present in the fixed point stateoriginating from different sources such as the compressor orelectronic noises. Therefore, the fixed point state with lowamplitude aperiodic noisy fluctuations would give H valuesbetween 0.2 and 0.5. During the transition to the state of limitcycle oscillations, H approaches zero. H captures the peri-odicity (or the loss of fractal nature) in the data even if theamplitudes are very low. We start observing an emergenceof periodicity as we approach the tipping point. The inher-ent noisy fluctuations perturb the system from the stable fixedpoint. These noise-induced oscillations occur at frequenciescentred around the natural frequency of the system and the os-cillations decay in time. These fluctuations contain very lowamplitude bursts of periodicity which has oscillations aroundthe natural frequency of the system . Capturing this slightperiodicity in the system variables close to the transition, H starts to decrease towards 0 well-before the rms of fluctua-tions grows. For very slow rates, H tends to decrease before µ (Fig. 4f). For relatively faster rates, the tipping is delayedsignificantly from µ ; nevertheless, H forewarns the tippingwell-before the rise in p (cid:48) rms (Fig. 4l,r).We repeat the computations of all EWS for data acquiredat many different rates of change of V . The results shown inFig. 4 are for three representative cases from this collection ofdata. A similar inference is obtained for even faster rates ofchange of voltage up to 240 mV/s. In all the cases, AC and H prove to be the better measures to forewarn an impending TAIcompared to other measures such as VAR, SKEW and K , andboth AC and H have comparable effectiveness.In summary, there are two things happening on approach-ing the transition; one is the growth of the amplitude of oscil-lations due to Hopf bifurcation, and the second is the increasein temporal correlation. Out of all the EWS discussed here, rms and variance capture the growth of amplitude of oscilla-tions; skewness and kurtosis detect the changing distributionof data during the transition; lag-1 autocorrelation determinesthe increasing correlation between consecutive time instants; H looks at the increasing correlation as well as the emergenceof periodicity. AC and H are capturing features that are dif-ferent from the features captured by p (cid:48) rms , VAR, SKEW and K . Variance, skewness and kurtosis do not change if we wereto shuffle the data randomly. However, AC and H will changeupon shuffling the data. Hence, AC and H are capturing thetemporal correlations present in the data more than just thestatistical characteristics of the data. From Fig. 5, it is clearthat the increase in correlation occurs prior to the growth ofamplitude. Therefore, AC and H are able to predict the tippingwell in advance. For faster rates of change of control param-eter, growth of correlation occurs at a much lower parametervalue than the growth of amplitude in the signal. C. Variation of warning time with rate
Our analysis shows that AC and H are effective EWS forcritical transitions in the considered thermoacoustic system.Next, we compare the early warning time for different ratesusing AC and H . Till now, we were focused on the warningin the parameter space, i.e. , at what parameter value we canpredict, compared to the tipping point. Ultimately, EWS needto be compared across both temporal domain and parameterspace.Figure 6-7 show the variation of H and AC as a function ofthe control parameter (left column) and time (right column),for different values of r . Variation of p (cid:48) rms is plotted on theright axis for comparison. A window of 20 s time interval isshown (right column) for all the plots, for the sake of com-parison. The time t (cid:48) is selected as a time instant before thetipping point within a window of 20 s. The warning time ismarked as τ (blue shaded region in Figs. 6-7). We calculate τ by selecting a threshold in AC and H so that, when it crossesthe threshold we are warned of an impending tipping. In thepresent case, we select AC = 0.4 and H = 0.1 as the thresh-old, because H never reduces to 0.1 unless the system is pro-ceeding towards an impending tipping. Similarly, AC also hasfluctuations, but, AC = 0.4 is out of bounds for the quiescentstate. The green dotted line in Figs. 6-7 corresponds to theselected threshold value of H = 0.1 or AC = 0.4. The time be-tween H = 0.1 or AC = 0.4 and the onset of TAI (correspondsto the maximum rate of change of p (cid:48) rms ) is considered as τ .This choice of threshold for AC and H are not unique, andany other slightly different threshold also will work equallywell. Before the amplitude rises steeply, we have more warn-ing time ( τ ) in the case of a relatively slower rates shown inFig. 6b,d.In terms of the parameter, we get warning at a parametervalue well ahead of the tipping point for faster rates (Figs. 6k,land 7k,l). However, the time needed to reach the tipping pointis really short, as the rate of change of P is fast. Hence, wehave relatively lesser time to implement control actions forfaster rates (Figs. 6k,l and 7k,l). For slower rates, we havemore time for control, but we are close to the tipping pointin the parameter space (Fig. 6a,b and 7a,b). As we are go-ing closer to the tipping point in case of slower rates, there isa possibility of N-tipping if external noise or inherent fluctu-ations of significant magnitude are present. The analysis ofwarning time using AC provides a similar inference, whereinslower rates result in more warning time and faster rates givemore early warning in terms of the control parameter (Fig. 7).These findings are summarized in Fig. 8.Next, we examine how the warning time ( τ ) varies with therate of change of voltage with time ( r = dV / dt ). We calcu-late τ by selecting a threshold, as mentioned in the precedingparagraph. As we increase the rate, τ decreases drastically(see Fig. 9). We observe an inverse power law relation be-tween the warning time and the rate of change of parameteras τ ∝ r n with n = - 0.76 ± n = -0.84 ± H , respectively. The uncertainty in fitting is estimatedwith 95% confidence. The scaling with a constant exponentholds for a range thresholds for AC (threshold: 0.2-0.8) and FIG. 6. Variation of p (cid:48) rms and H evaluated for data acquired at different rates of change of V is plotted as a function of the P (left column -a,c,e,g,i,k) and time (right column - b,d,f,h,j,l). Rate of change of parameter increases from top to bottom. For all these rates, we are able topredict the tipping using H , before p (cid:48) rms starts to increase. The green dotted line represents to the threshold value of H = 0.1. The warningtime before the tipping (marked with blue colour) is the difference between the time at which H crosses the threshold and the time at whichmaximum rate of change of p (cid:48) rms is observed. H (threshold: 0.15-0.06). The scaling using different thresh-olds are shown in Appendix B. D. Model for noisy Hopf bifurcation with continuousvariation of parameter
In the present study, we use the model of a nonlin-ear oscillator with additive noise exhibiting subcritical Hopfbifurcation .¨ η + α ˙ η + ω η = ˙ η (cid:0) β + K η − γη (cid:1) + N , (1) d α dt = rt , where η and ˙ η are the state variables, α and β are lineardamping and driving respectively and ω is the natural fre-quency. We use additive Gaussian white noise N of intensity Γ with an autocorrelation < N N τ > = Γ δ ( τ ) . FollowingNoiray , the values of the parameters ω , β , γ and K are keptconstant as follows: ω = 2 π ×
120 rad/s, β = 50 rad/s, γ = 0.7,and K = 9. We select a low value of Γ = √ − as we aremodelling a laminar system. Here, the linear damping ( α ) canbe considered as the bifurcation parameter analogous to thecontrol parameter in the experiment. α is reduced from 200rad/s to -50 rad/s to capture the subcritical Hopf bifurcation FIG. 7. Variation of p (cid:48) rms and AC is plotted as a function of control parameter, P , (left column - a,c,e,g,i,k) and time, t , (right column -b,d,f,h,j,l) evaluated for data acquired at different rates of change of P . Rate of change of control parameter increases from top to bottom. Forall the transitions shown, AC provides an early warning before p (cid:48) rms starts to rise. The green dotted line represents to the threshold value of AC .The warning time, τ , (marked with blue colour) is the difference between the time at which AC crosses the threshold and the time at whichmaximum rate of change of p (cid:48) rms is observed. from a state of stable fixed point to limit cycle oscillations.In order to study the effect of rate of change of parameterin the model, we vary α continuously. In experiments, theheater power is changing nonlinearly even though we vary V linearly. To replicate a similar variation of the control param-eter, we change α as α = α + r t in the model. Initially, thelinear damping α is high (200 rad/s) to obtain a fixed point.Upon decreasing α , once α = β , the fixed point becomes un-stable. However, there has to be a non-zero perturbation forthe system to escape from the fixed point. The additive whitenoise term ( N ) in Eq. 1 constantly perturbs the system fromthe fixed point and helps to jump to the limit cycle state. In ex-periments, even though the system is laminar, low amplitudenoisy fluctuations are always present in the flow field.Although the above model captures the dynamic transi-tions qualitatively, the initial fixed point state appears to havebursts with a high level of periodicity (Fig. 10b), unlike the noisy aperiodic fluctuations observed in the experiments. No-tably, the fixed point in Hopf bifurcation under the influ-ence of noise contains fluctuations with bursts of periodic-ity with shallow peaks near the fundamental frequency of thesystem . However, the experimental data appears to bemore aperiodic. In experiments, the base state with constantairflow and zero heater power itself generates low amplitudeaperiodic pressure fluctuations which will also get measuredalong with the system dynamics. This inherent fluctuationsand the measurement noise involved in the dynamics could bemodelled by adding Gaussian white noise ( N ) with intensity Γ to the time series of η ( t ) obtained by solving the model. η ( t ) = η ( t ) + N ( t ) . (2)We choose a noise intensity ratio, Γ / Γ = 0.03, to repli-cate the experimental data qualitatively. A very low value of Γ / Γ would result in high values of AC and low values of H FIG. 8. The effect of rate of change of control parameter on EWS for critical transitionsFIG. 9. Warning time ( τ ) obtained from (a) AC and (b) H are plottedas a function of r in logarithmic scale. (b) The value of τ decreaseswith r following an inverse power law relation as τ ∝ r n with n =-0.76 ± n = -0.84 ± H , respectively. Thethreshold for obtaining warning is 0.4 for AC and 0.1 for H . due to the presence of low amplitude bursty oscillations dur-ing the noisy fixed point state. As we gradually increase theintensity of external noise ( i.e., increasing Γ / Γ ), the valuesof H increases and AC decreases. Then, the values during thetransition also matches with those obtained from experiments.However, if we add external noise more than necessary, itwould suppress all the underlying dynamics and produce AC= 0 and H = . α is shown in Fig. 11. For all the measures, we observe a sim-ilar performance as that observed for the experiments. Again, FIG. 10. (a) Bifurcation diagram showing the variation of η (cid:48) rms as afunction of linear damping, α . The Hopf point ( µ ) is at α = 50 rad/s,where α = β . (b) The time series of η obtained by solving Eq. 1 fora continuous variation of α is shown along with a zoomed view ofthe bursts occurring during the fixed point state. AC and H provide a better warning for an impending TAIcompared to the other measures. For faster rates of changeof parameter, SKEW and K predict the tipping slightly betterfor the model when compared to the experiments. For slowerrates of change of parameter, the transition is more abrupt interms of the control parameter and all measures except ACand H detect this very close to the tipping point. In contrast,the growth of amplitude of fluctuations is more gradual forfaster rates of change of parameter, and the other measures(SKEW and K ) are also able to predict the tipping before η (cid:48) rms and VAR increases.Next, we analyze the warning time obtained using AC and H for different rates of change of α . The warning time re-duces with the rate of change of parameter following an in-1 FIG. 11. Variation of η (cid:48) rms , lag-1 autocorrelation (AC), variance (VAR), skewness (SKEW), kurtosis ( K ) and Hurst exponent ( H ) with theparameter α during the transition to limit cycle oscillations in the model. Each column corresponds to the results for a particular rate of changeof α ((a)-(f): r = -0.0005, (g)-(l): r = -0.05, (m)-(r): - r = -5). The point of the maximum rate of change of η (cid:48) rms is marked with a red dottedline. For higher values of r , kurtosis and skewness detect the tipping before η (cid:48) rms and VAR. Here, AC and H consistently give early warningwell-before the transition to limit cycle oscillations for all the rates.FIG. 12. Warning time ( τ ) obtained from (a) AC and (b) H are plotted as a function of r in a double logarithmic plot for both linear (bluecolour) and quadratic (green colour) variation of α . We observe a power law relation between τ and AC and τ and H . The threshold chosenfor obtaining warning is 0.4 for AC and 0.1 for H . α = α + r l t and α = α + r nl t (Fig. 12). The trends are qualitatively thesame, even though the obtained exponents are different fromthose of experiments. This difference in the exponent couldbe because of multiple reasons. The transition depends on therate of change of parameter, the initial value of the parame-ter and the type and intensity of noise present in the system.Also, the other parameters such as the β , κ and γ can be var-ied to adjust the model to represent the experimental systembetter. In the current study, we vary V linearly; ideally, P hasto vary as P = V ( t ) / R , where R is the resistance. However, R can change slightly with temperature as V increases. In realsystems, as we vary one parameter in the system, there canbe multiple parameters changing simultaneously. Hence, theactual control parameter in practical systems can be a combi-nation of multiple parameters. We do not attempt to obtain theexact parameters and condition, but only aim for a qualitativematch between the experiment and the model. E. Relation between lag-1 autocorrelation and variance closeto the tipping point
Critical slowing down leads to an increased autocorrelationand variance of fluctuations in a stochastically perturbed sys-tem approaching a bifurcation . Near the critical point, ACtends to one, and VAR tends to infinity. Here, we examine thevariation of lag-1 autocorrelation (AC) with the variance offluctuations close to the tipping where the AC grows and sat-urates to 1. VAR stays nearly zero almost till the onset of TAIand then shoots up to really high values. At the same time,AC increases gradually during the transition. Here, we ob-serve that AC scales with VAR following a hyperexponentialrelation , d ( VAR ) d ( AC ) = k ( VAR ) a . (3)Refer to Fig. 13a for the scaling obtained from experimentsand Fig. 13b for data from the model. The most commongrowth principle is exponential, where a = a , we fit d ( VAR ) d ( AC ) with VAR in a double logarithmic plot (Figs. 13c-d). The fitting isperformed only for the data points in the middle of the curvein Fig. 13a. The scaling shown in Fig. 13c is a representa-tive case of one experiment with a rate of 3 mV/s. The expo-nent is found to be a = 1.94 ± a = 1.99 ± VAR = [ a ( k AC + constant )] ( / a ) , where a = − a + a ∼ ∝ − k AC + constant . If we apply the limitsapproaching the tipping point: VAR tends to infinity as ACtends to one, we get the constant as − k and the final expres-sion as, VAR ∝ − AC . (4) We observe the same scaling exponent for linear and nonlin- FIG. 13. (a-b) Increase in the variance of fluctuations as a functionof AC shows a curve increasing faster than exponential (hyperexpo-nential). (c-d) To find the exponent of the hyperexponential function,a linear fit is done on the double logarithmic plot of d ( A ) d ( AC ) vs . VAR. ear variation of the control parameter in the model. The scal-ing relation seems to be independent of the functional form ofthe parameter. This hyperexponential scaling observed in ex-periments and model irrespective of the rate of change of pa-rameter and the functional form of control parameter suggeststhat this can be a scaling during the occurrence of dynamicHopf bifurcation. IV. CONCLUSION
We study several early warning signals for critical transi-tions in a thermoacoustic system. Compared to the quasi-static bifurcation, the onset of tipping is delayed when thecontrol parameter is varied continuously at a finite rate. Weconfirm the observation of increased delay with increase inthe rate of change of control parameter. By analyzing the per-formance of various early warning signals, we observe thatthe variance, kurtosis and skewness do not provide adequatewarning; they change only when p (cid:48) rms rises. The lag-1 au-tocorrelation and the Hurst exponent are able to predict thetransition well-before the tipping point. We confirmed thisobservation by performing experiments at different rates ofchange of control parameter. For slower rates, AC and H givemore warning time compared to faster rates, even though weare relatively close to the transition in terms of the parameter.On the other hand, for faster rates where we have relativelylesser time to initiate control measures, AC and H capturethe tipping at a parameter value which is well ahead of thetipping point. Furthermore, we notice that the warning timereduces with the rate of change of parameter following an in-verse power law relation. Then, we perform a similar analysisfor a noisy Hopf bifurcation model. The qualitative featuresof the EWS using the lag-1 autocorrelation and the Hurst ex-3ponent are captured using the model. Finally, we empiricallyobtained a relation between lag-1 autocorrelation and the vari-ance of fluctuations for dynamic Hopf bifurcation. This hy-perexponential scaling is found to be independent of the func-tional form of variation of the control parameter . ACKNOWLEDGMENTS
We acknowledge the discussions and help from Dr. VishnuR. Unni, Dr. Gopalakrishnan and Dr. Samadhan Pawar. I. P.is grateful to the Ministry of Human Resource Development,India and Indian Institute of Technology Madras for providingresearch assistantship. R. I. S. acknowledges the funding fromJ.C. Bose fellowship (JCB/2018/000034/SSC).
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request.
Appendix A: Methodology of computation of early warningmeasures
We provide an overview of the different EWS used in thepresent study.
1. Root mean square
The root mean square ( rms ) of a time series, p ( t ) , is definedas the square root of the mean square or the arithmetic meanof the squares of all the elements in the time series. It is alsoknown as the quadratic mean. p rms = (cid:118)(cid:117)(cid:117)(cid:116) N (cid:32) N ∑ i = p i (cid:33) (A1)where N is the total number of data points. Root mean squarevalue has been widely used in many engineering applicationsas a first step to check the sudden increase in the amplitude offluctuations about the mean.
2. EWS based on Critical slowing down
Systems approaching a transition where the current statebecomes unstable and transitions to another state shows slowresponse to small external perturbations. This phenomenonof slow return rate is known as critical slowing down andcan be detected by increased autocorrelation and variance offluctuations . a. Autocorrelation Autocorrelation is the correlation of a signal, p ( t ) , with adelayed copy of itself as a function of delay ( τ ). It is definedas follows: AC ( τ ) = ∑ Ni = p ( i ) p ( i − τ ) σ (A2)In the current study, we consider lag-1 autocorrelation (AC)which computes the correlation between values that are onetime step apart. Critical slowing down leads to an increasein the short term memory of a system and can captured bycorrelation at low lags. b. Variance Variance (VAR) is the expectation of the squared differ-ences from its mean. VAR measures how the data is spreadout from their average value and it is the second moment ofthe distribution. VAR = N N ∑ i = ( p i − p ) (A3)where p is the mean and N is the number of data points.
3. Skewness and kurtosis
Skewness (SKEW) is a measure of symmetry or the lackof symmetry of the probability distribution of data about itsmean. A distribution is symmetric if it looks the same to theleft and right of the mean value. The skewness for a normaldistribution is zero and skewness will be near zero for anysymmetric data. The skewness of a random variable X is itsthird moment. SKEW = ∑ Ni = ( p i − p ) / N σ (A4)where p is the mean, σ is the standard deviation, and N isthe number of data points. Negative SKEW indicates that thedistribution is skewed left and positive SKEW indicates thatthe distribution is skewed right. Skewed left means that theleft tail is longer than the right tail and vice versa. Close to atipping point, probability distributions may become asymmet-ric with a non-zero skewness . Kurtosis ( K ) gives informa-tion about whether the probability distribution of the data hasheavy tails, or is more centred compared to a normal distri-bution. The kurtosis of a normal distribution is 3, and higherkurtosis indicates more outliers or heavy tail in the data. IfKurtosis is less than 3, it means that the distribution has lesseroutliers compared to the normal distribution. The kurtosis isthe fourth standardized moment and is defined as follows: K = ∑ Ni = ( p i − p ) / N σ (A5)where p is the mean, σ is the standard deviation and N is thenumber of data points.4 FIG. 14. The inverse power law scaling between the warning time and rate of change of parameter is plotted for different values of thresholdsfor (a) AC and (b) H .
4. Hurst exponent
The Hurst exponent ( H ) is related to the fractal dimension( D ) as D = − H , in the case of a time series. Generally, H hasvalues between 0 and 1 for a time series of fractal dimensionbetween 1 and 2. We calculate H following the algorithm ofMultifractal detrended fluctuation analysis . Multifractal detrended fluctuation analysis (MFDFA)
First, we take a time series x ( t ) of length N . The mean sub-tracted cumulative deviate series Y ( k ) can be defined as, Y ( k ) = k ∑ t = [ x t − ¯ x ] , k = , , ..., N , (A6)where ¯ x is the mean of the time series. We divide the devi-ate series Y ( k ) into N w = [ N / w ] non-overlapping segments ofequal length w , where [ N / w ] represents the greatest integerfunction. Then, we find a polynomial fit ( ¯ Y i ) for each of thesesegments i and we subtract the polynomial fit from the deviateseries ( Y i ) to obtain the fluctuations. In the present study, weuse a polynomial fit of order 1. Then, the variance of fluctua-tions is determined as, F ( w , i ) = w (cid:104) w ∑ t = ( Y i ( t ) − ¯ Y i ) (cid:105) , (A7)for each segment i = , , ..., N w .The structure function of order 2 and span w , F w is defined asfollows: F w = (cid:104) N w N w ∑ i = F ( w , i ) (cid:105) / . (A8)We repeat the same steps for different time scales or span w and plot the variation of F w with the span w in a logarithmicscale. The slope of the linear regime for a range of span sizes w is known as the Hurst exponent ( H ). Appendix B: Robustness of EWS with the threshold
We observe that the inverse power law scaling between thewarning time and rate of change of parameter is consistentwith almost the same exponents for different values of thresh-olds of EWS as shown in Fig. 14.
FIG. 15. Acoustic pressure fluctuations ( p (cid:48) ) acquired for a constantvalue of control parameter (voltage = 1.6 V and heater power = 253.7W), lower than the Hopf point. The time series contains only lowamplitude aperiodic fluctuations, and there is no transition to limitcycle oscillations. The rms of pressure fluctuations is plotted alongwith the signal in the top plot. The corresponding variation of allthe EWS is shown as subplots. We calculate the EWS for a movingwindow of 1 s with an overlap of 0.98 s. p (cid:48) rms , AC, VAR, SKEW, K and H stay nearly constant for the total duration of the experiment,and we do not observe any significant change in the values of EWS. Appendix C: Analysis to check false warnings
To check for false warnings, we calculate EWS for data forcases where transition to thermoacoustic instability does notoccur. Here, we have the time series data acquired for con-stant values of voltage (quasi-static experiments). Figure 15shows one such data and the corresponding variation of all theEWS. This is a representative case of data for voltage = 1.6V and heater power = 253.7 W, and we have confirmed theseobservations by analysing many data sets which are taken forquasi-static experiments. We wait at a particular value of con-trol parameter far below the Hopf point, and there is no tran-sition to limit cycle oscillations as expected. All the EWScalculated for a moving window show constant values indica-tive of a low amplitude aperiodic state, for the entire lengthof the data. Hence, we do not observe any false warnings forthese EWS. T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf,and H. J. Schellnhuber, “Tipping elements in the earth’s climate system,”Proceedings of the national Academy of Sciences , 1786–1793 (2008). M. Scheffer, S. Carpenter, J. A. Foley, C. Folke, and B. Walker, “Catas-trophic shifts in ecosystems,” Nature , 591 (2001). D. Sornette and A. Johansen, “Large financial crashes,” Physica A: Statis-tical Mechanics and its Applications , 411–422 (1997). J. G. Venegas, T. Winkler, G. Musch, M. F. V. Melo, D. Layfield,N. Tgavalekos, A. J. Fischman, R. J. Callahan, G. Bellani, and R. S. Har-ris, “Self-organized patchiness in asthma as a prelude to catastrophic shifts,”Nature , 777 (2005). B. Litt, R. Esteller, J. Echauz, M. D’Alessandro, R. Shor, T. Henry, P. Pen-nell, C. Epstein, R. Bakay, M. Dichter, et al. , “Epileptic seizures may beginhours in advance of clinical onset: a report of five patients,” Neuron ,51–64 (2001). P. E. McSharry, L. A. Smith, and L. Tarassenko, “Prediction of epilepticseizures: are nonlinear methods relevant?” Nature medicine , 241 (2003). S. R. Carpenter, D. Ludwig, and W. A. Brock, “Management of eutroph-ication for lakes subject to potentially irreversible change,” Ecological ap-plications , 751–771 (1999). National-Research-Councils and other,
New directions for understandingsystemic risk: a report on a conference cosponsored by the Federal Re-serve Bank of New York and the National Academy of Sciences (NationalAcademies Press, 2007). R. M. May, S. A. Levin, and G. Sugihara, “Ecology for bankers,” Nature , 893–894 (2008). P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, “Tipping points in opensystems: bifurcation, noise-induced and rate-dependent examples in the cli-mate system,” Philosophical Transactions of the Royal Society A: Mathe-matical, Physical and Engineering Sciences , 1166–1184 (2012). J. Tony, S. Subarna, K. S. Syamkumar, G. Sudha, S. Akshay, E. A.Gopalakrishnan, E. Surovyatkina, and R. I. Sujith, “Experimental investi-gation on preconditioned rate induced tipping in a thermoacoustic system,”Scientific reports , 1–7 (2017). S. M. Baer, T. Erneux, and J. Rinzel, “The slow passage through a hopfbifurcation: delay, memory effects, and resonance,” SIAM Journal on Ap-plied mathematics , 55–71 (1989). G. Bonciolini, D. Ebi, E. Boujo, and N. Noiray, “Experiments and mod-elling of rate-dependent transition delay in a stochastic subcritical bifurca-tion,” Royal Society open science , 172078 (2018). Y. Park, Y. Do, and J. M. Lopez, “Slow passage through resonance,” Phys-ical Review E , 056604 (2011). N. Berglund, “Dynamic bifurcations: Hysteresis, scaling laws and feedbackcontrol,” Progress of Theoretical Physics Supplement , 325–336 (2000). A. Majumdar, J. Ockendon, P. Howell, and E. Surovyatkina, “Transitionsthrough critical temperatures in nematic liquid crystals,” Physical ReviewE , 022501 (2013). V. R. Unni, E. A. Gopalakrishnan, K. S. Syamkumar, R. I. Sujith,E. Surovyatkina, and J. Kurths, “Interplay between random fluctuations and rate dependent phenomena at slow passage to limit-cycle oscillationsin a bistable thermoacoustic system,” Chaos: An Interdisciplinary Journalof Nonlinear Science , 031102 (2019). M. Schroeder,
Fractals, chaos, power laws: Minutes from an infinite par-adise (Courier Corporation, 2009). M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter,V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara, “Early-warning signals for critical transitions,” Nature , 53 (2009). V. Dakos, M. Scheffer, E. H. van Nes, V. Brovkin, V. Petoukhov, andH. Held, “Slowing down as an early warning signal for abrupt climatechange,” Proceedings of the National Academy of Sciences , 14308–14312 (2008). T. Wilkat, T. Rings, and K. Lehnertz, “No evidence for critical slowingdown prior to human epileptic seizures,” Chaos: An Interdisciplinary Jour-nal of Nonlinear Science , 091104 (2019). T. T. Tsotsis, R. C. Sane, and T. H. Lindstrom, “Bifurcation behavior ofa catalytic reaction due to a slowly varying parameter,” AIChE journal ,383–388 (1988). A. K. Kapila, “Arrhenius systems: dynamics of jump due to slow pas-sage through criticality,” SIAM Journal on Applied Mathematics , 29–42(1981). L. M. Bilinsky and S. M. Baer, “Slow passage through a hopf bifurcation inexcitable nerve cables: Spatial delays and spatial memory effects,” Bulletinof mathematical biology , 130–150 (2018). P. Ashwin, C. Perryman, and S. Wieczorek, “Parameter shifts for nonau-tonomous systems in low dimension: bifurcation-and rate-induced tipping,”Nonlinearity , 2185 (2017). W. Scharpf, M. Squicciarini, D. Bromley, C. Green, J. Tredicce, and L. Nar-ducci, “Experimental observation of a delayed bifurcation at the thresholdof an argon laser,” Optics communications , 344–348 (1987). A. N. Pisarchik, R. Jaimes-Reátegui, C. A. Magallón-García, and C. O.Castillo-Morales, “Critical slowing down and noise-induced intermittencyin bistable perception: bifurcation analysis,” Biological Cybernetics ,397–404 (2014). K. I. Matveev,
Thermoacoustic instabilities in the Rijke tube: experimentsand modeling , Ph.D. thesis, California Institute of Technology (2003). M. P. Juniper, “Triggering in the horizontal rijke tube: non-normality, tran-sient growth and bypass transition,” Journal of Fluid Mechanics , 272–308 (2011). E. A. Gopalakrishnan and R. I. Sujith, “Effect of external noise on the hys-teresis characteristics of a thermoacoustic system,” Journal of Fluid Me-chanics , 334–353 (2015). M. P. Juniper and R. I. Sujith, “Sensitivity and nonlinearity of thermoacous-tic oscillations,” Annual Review of Fluid Mechanics , 661–689 (2018). S. C. Fisher and S. A. Rahman, “Remembering the giants: Apollo rocketpropulsion development,” (2009). T. C. Lieuwen and V. Yang,
Combustion instabilities in gas turbine engines:operational experience, fundamental mechanisms, and modeling (Ameri-can Institute of Aeronautics and Astronautics, 2005). C. Fleming, “Turbine makers are caught in innovation trap,” Wall StreetJournal (February 13, 1998). S. H. Strogatz, M. Friedman, A. J. Mallinckrodt, S. McKay, et al. , “Nonlin-ear dynamics and chaos: With applications to physics, biology, chemistry,and engineering,” Computers in Physics , 532–532 (1994). V. Dakos, E. H. Van Nes, P. D’Odorico, and M. Scheffer, “Robustness ofvariance and autocorrelation as indicators of critical slowing down,” Ecol-ogy , 264–271 (2012). T. M. Lenton, R. J. Myerscough, R. Marsh, V. N. Livina, A. R. Price, S. J.Cox, and GENIE-team, “Using GENIE to study a tipping point in the cli-mate system,” Philosophical Transactions of the Royal Society A: Mathe-matical, Physical and Engineering Sciences , 871–884 (2009). V. Guttal and C. Jayaprakash, “Changing skewness: an early warning signalof regime shifts in ecosystems,” Ecology letters , 450–460 (2008). V. Nair and R. I. Sujith, “Multifractality in combustion noise: predicting animpending combustion instability,” Journal of Fluid Mechanics , 635–655 (2014). R. I. Sujith and V. R. Unni, “Complex system approach to investigate andmitigate thermoacoustic instability in turbulent combustors,” Physics ofFluids , 061401 (2020). B. B. Mandelbrot,
The fractal geometry of nature , Vol. 173 (WH freemanNew York, 1983). H. E. Hurst, “Long-term storage capacity of reservoirs,” Trans. Amer. Soc.Civil Eng. , 770–799 (1951). J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin,A. Bunde, and H. E. Stanley, “Multifractal detrended fluctuation analy-sis of nonstationary time series,” Physica A: Statistical Mechanics and itsApplications , 87–114 (2002). P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosen-blum, Z. R. Struzik, and H. E. Stanley, “Multifractality in human heartbeatdynamics,” Nature , 461 (1999). K. Hu, P. C. Ivanov, M. F. Hilton, Z. Chen, R. T. Ayers, H. E. Stanley, andS. A. Shea, “Endogenous circadian rhythm in an index of cardiac vulner-ability independent of changes in behavior,” Proceedings of the NationalAcademy of Sciences , 18223–18227 (2004). D. Grech and Z. Mazur, “Can one make any crash prediction in financeusing the local Hurst exponent idea?” Physica A: Statistical Mechanics andits Applications , 133–145 (2004). N. Vandewalle and M. Ausloos, “Coherent and random sequences in fi-nancial fluctuations,” Physica A: Statistical Mechanics and its Applications , 454–459 (1997). D. Grech and G. Pamuła, “The local Hurst exponent of the financial time se-ries in the vicinity of crashes on the polish stock exchange market,” PhysicaA: statistical mechanics and its applications , 4299–4308 (2008). J. Alvarez-Ramirez, J. Alvarez, E. Rodriguez, and G. Fernandez-Anaya,“Time-varying Hurst exponent for us stock markets,” Physica A: StatisticalMechanics and its Applications , 6159–6169 (2008). J. A. Matos, S. M. Gama, H. J. Ruskin, A. Al Sharkasi, and M. Crane,“Time and scale Hurst exponent analysis for financial markets,” Physica A:Statistical Mechanics and its Applications , 3910–3915 (2008). K. Domino, “The use of the Hurst exponent to predict changes in trendson the warsaw stock exchange,” Physica A: Statistical Mechanics and itsApplications , 98–109 (2011). V. Suyal, A. Prasad, and H. P. Singh, “Nonlinear time series analysis ofsunspot data,” Solar Physics , 441–449 (2009). A. Kilcik, C. N. K. Anderson, J. P. Rozelot, H. Ye, G. Sugihara, andA. Ozguc, “Nonlinear prediction of solar cycle 24,” The Astrophysical Jour-nal , 1173 (2009). V. R. Unni and R. I. Sujith, “Multifractal characteristics of combustor dy-namics close to lean blowout,” Journal of Fluid Mechanics , 30–50(2015). H. Gotoda, M. Amano, T. Miyano, T. Ikawa, K. Maki, and S. Tachibana,“Characterization of complexities in combustion instability in a lean pre-mixed gas-turbine model combustor,” Chaos: An Interdisciplinary Journalof Nonlinear Science , 043128 (2012). S. Havlin, L. A. N. Amaral, Y. Ashkenazy, A. L. Goldberger, P. C. Ivanov,C.-K. Peng, and H. E. Stanley, “Application of statistical physics to heart-beat diagnosis,” Physica A: Statistical Mechanics and its Applications ,99–110 (1999). L. Telesca, V. Cuomo, V. Lapenna, and M. Macchiato, “A new approachto investigate the correlation between geoelectrical time fluctuations andearthquakes in a seismic area of southern italy,” Geophysical Research Let-ters , 4375–4378 (2001). B. Qian and K. Rasheed, “Hurst exponent and financial market predictabil-ity,” in
IASTED conference on Financial Engineering and Applications (2004) pp. 203–209. E. A. Gopalakrishnan, Y. Sharma, T. John, P. S. Dutta, and R. I. Sujith,“Early warning signals for critical transitions in a thermoacoustic system,”Scientific reports , 1–10 (2016). G. Ghanavati, P. D. H. Hines, T. I. Lakoba, and E. Cotilla-Sanchez, “Un-derstanding early indicators of critical transitions in power systems fromautocorrelation functions,” IEEE Transactions on Circuits and Systems I:Regular Papers , 2747–2760 (2014). K. Wiesenfeld, “Noisy precursors of nonlinear instabilities,” Journal of Sta-tistical Physics , 1071–1097 (1985). N. Noiray, “Linear growth rate estimation from dynamics and statistics ofacoustic signal envelope in turbulent combustors,” Journal of Engineeringfor Gas Turbines and Power , 041503 (2017). H. Fujisaka and M. Inoue, “Correlation functions of temporal fluctuations,”Physical Review A , 4778 (1989). S. D. Varfolomeyev and K. G. Gurevich, “The hyperexponential growth ofthe human population on a macrohistorical scale,” Journal of TheoreticalBiology , 367–372 (2001). E. A. F. E. Ihlen, “Introduction to multifractal detrended fluctuation analysisin Matlab,” Frontiers in physiology3