Science, Faculty of
Mathematics, Department of
DSpace
UBCV
Zhu, Jielin
Doctor of Philosophy - PhD
University of British Columbia
We consider the effect on tipping from an additive periodic forcing and an additive white noise in a canonical model with a saddle-node bifurcation and a slowly varying bifurcation parameter. Using a multiple scales analysis, we consider the effect of amplitude and frequency of the periodic forcing relative to the drifting rate of the varying bifurcation parameter. We show that a high frequency oscillation drives an earlier tipping when the bifurcation parameter varies more slowly, with the advance of the tipping point proportional to the square of the ratio of amplitude to frequency. In the low frequency case the position of the tipping point is affected by the frequency, amplitude and phase of the low frequency oscillation. The results are based on an analysis of the local concavity of the trajectory. The tipping point location increases with the amplitude of the periodic forcing, with critical amplitudes where there are jumps in the location. Using a WKB method and a multi-scale analysis, we consider the effect of the amplitude of the noise relative to the drifting rate of the bifurcation parameter. We show that the early tipping is likely to happen if the drifting rate of the bifurcation parameter is smaller than the square of the amplitude of the noise. The WKB approximation shows that the quasi-stationary PDF is related to the potential well of the system. The probability of early tipping, which is related to the PDF near the unstable equilibrium, should be considered locally through a convection-diffusion process. The analysis of the Morris-Lecar model with oscillations shows that the method of multiple scales we used is applicable to the higher dimensional system. For the sea ice model, the method of multiple scales is adapted for the case when the amplitude of the oscillation is larger than its frequency. We applied the WKB method and the multi-scale analysis on the zero-dimensional climate model which has an asymmetric structure of the saddle-node bifurcation. The approximation of the PDF shows that the asymmetry of the shape of the bifurcation diagram plays an important role in the probability of early tipping.
ASYMPTOTIC BEHAVIOUR OFTIPPING POINTS INNON-AUTONOMOUS SYSTEMSWITH RANDOM AND PERIODICFORCINGbyJielin ZhuB.Sc., Shandong University, 2007M.Sc., Trent University, 2009A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2015c© Jielin Zhu 2015AbstractWe consider the effect on tipping from an additive periodic forcing and anadditive white noise in a canonical model with a saddle-node bifurcationand a slowly varying bifurcation parameter. Using a multiple scales analy-sis, we consider the effect of amplitude and frequency of the periodic forcingrelative to the drifting rate of the varying bifurcation parameter. We showthat a high frequency oscillation drives an earlier tipping when the bifur-cation parameter varies more slowly, with the advance of the tipping pointproportional to the square of the ratio of amplitude to frequency. In the lowfrequency case the position of the tipping point is affected by the frequency,amplitude and phase of the low frequency oscillation. The results are basedon an analysis of the local concavity of the trajectory. The tipping pointlocation increases with the amplitude of the periodic forcing, with criticalamplitudes where there are jumps in the location. Using a WKB methodand a multi-scale analysis, we consider the effect of the amplitude of thenoise relative to the drifting rate of the bifurcation parameter. We showthat the early tipping is likely to happen if the drifting rate of the bifur-cation parameter is smaller than the square of the amplitude of the noise.The WKB approximation shows that the quasi-stationary PDF is relatedto the potential well of the system. The probability of early tipping, whichiiAbstractis related to the PDF near the unstable equilibrium, should be consideredlocally through a convection-diffusion process. The analysis of the Morris-Lecar model with oscillations shows that the method of multiple scales weused is applicable to the higher dimensional system. For the sea ice model,the method of multiple scales is adapted for the case when the amplitude ofthe oscillation is larger than its frequency. We applied the WKB method andthe multi-scale analysis on the zero-dimensional climate model which has anasymmetric structure of the saddle-node bifurcation. The approximation ofthe PDF shows that the asymmetry of the shape of the bifurcation diagramplays an important role in the probability of early tipping.iiiPrefaceAll of the work presented in this dissertation is conducted by the authorJielin Zhu under the guidance of the supervisor Rachel Kuske.Chapter 2, Section 4.1 and Section 4.2 are based on work submitted [20].Zhu J., R. Kuske and T. Erneux, Tipping points near a delayed saddle-nodebifurcation with periodic forcing, SIAM J on Applied Dynamical Systems,Submitted. I conducted most of the analytical results and wrote most ofthe manuscript.Chapter 3 and Section 4.3 are original, unpublished work.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . 51.3 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . 102 Periodic Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 High Frequency Oscillation Ω ≫ 1 . . . . . . . . . . . . . . . 16vTable of Contents2.2.1 Constant bifurcation parameter . . . . . . . . . . . . 182.2.2 Slowly varying bifurcation parameter . . . . . . . . . 212.2.3 Large amplitude (A≫ 1) . . . . . . . . . . . . . . . . 292.3 Low Frequency Oscillation Ω ≪ 1 . . . . . . . . . . . . . . . 312.3.1 Ω = O(µ) or smaller . . . . . . . . . . . . . . . . . . . 312.3.2 Ω = µν (0 < ν < 1) . . . . . . . . . . . . . . . . . . . 413 Random Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . 473.1 Scaling Relationship . . . . . . . . . . . . . . . . . . . . . . . 473.2 WKB Approximation . . . . . . . . . . . . . . . . . . . . . . 513.3 PDF Tail Approximation . . . . . . . . . . . . . . . . . . . . 594 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.1 Morris-Lecar Model . . . . . . . . . . . . . . . . . . . . . . . 684.1.1 High frequency oscillation Ω ≫ 1 . . . . . . . . . . . 734.1.2 Low frequency oscillation Ω ≪ 1 . . . . . . . . . . . . 774.2 Energy Balance Sea Ice Model . . . . . . . . . . . . . . . . . 824.3 The Zero-dimensional Climate Model . . . . . . . . . . . . . 954.3.1 The model . . . . . . . . . . . . . . . . . . . . . . . . 954.3.2 Scaling relationship (asymmetry) . . . . . . . . . . . 984.3.3 Tail approximation (asymmetry) . . . . . . . . . . . . 1044.3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 1065 Conclusion and Future Work . . . . . . . . . . . . . . . . . . 119Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127viTable of ContentsAppendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .A Local Analysis for High Frequency Forcing . . . . . . . . . 133B Local Analysis for Low Frequency Forcing . . . . . . . . . . 136C Local Analysis for Amplitude Threshold . . . . . . . . . . . 138D Local Analysis for ML Model . . . . . . . . . . . . . . . . . . 141E Gaussian Approximation . . . . . . . . . . . . . . . . . . . . . 143viiList of Tables3.1 Summary of Scaling . . . . . . . . . . . . . . . . . . . . . . . 523.2 Summary of Variables . . . . . . . . . . . . . . . . . . . . . . 57viiiList of Figures2.1 Delay Effect of the Drifting . . . . . . . . . . . . . . . . . . . 172.2 Bifurcation Structure . . . . . . . . . . . . . . . . . . . . . . . 192.3 Periodic Solutions for Low Frequency . . . . . . . . . . . . . . 222.4 High-frequency Example . . . . . . . . . . . . . . . . . . . . . 262.5 High-frequency Tipping . . . . . . . . . . . . . . . . . . . . . 282.6 Low-frequency Example . . . . . . . . . . . . . . . . . . . . . 332.7 Low-frequency Threshold . . . . . . . . . . . . . . . . . . . . 362.8 Low-frequency Concavity . . . . . . . . . . . . . . . . . . . . 392.9 Low-frequency Amplitude Threshold . . . . . . . . . . . . . . 402.10 Low-frequency General Tipping . . . . . . . . . . . . . . . . . 442.11 Periodic Forcing Summary . . . . . . . . . . . . . . . . . . . . 453.1 Potential Well . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.2 PDF Comparison . . . . . . . . . . . . . . . . . . . . . . . . . 563.3 Merged PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.4 Modified PDF Comparison . . . . . . . . . . . . . . . . . . . 664.1 ML Bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . 714.2 ML High Frequency Case . . . . . . . . . . . . . . . . . . . . 76ixList of Figures4.3 ML Low Frequency Case . . . . . . . . . . . . . . . . . . . . . 804.4 Sea Ice Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.5 Sea Ice Tipping . . . . . . . . . . . . . . . . . . . . . . . . . . 914.6 Sea Ice Tipping Changes . . . . . . . . . . . . . . . . . . . . . 934.7 Asymmetric Tipping . . . . . . . . . . . . . . . . . . . . . . . 994.8 Quasi-stationary Solution . . . . . . . . . . . . . . . . . . . . 1074.9 Full Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 1084.10 Probability Comparison . . . . . . . . . . . . . . . . . . . . . 1104.11 Probability Comparison with Different Parameter Values . . . 1134.12 Skewness Comparison . . . . . . . . . . . . . . . . . . . . . . 114xAcknowledgementsI would like to record my gratitude and appreciation to my supervior, Dr.Rachel Kuske. Her guidance and advice from the beginning are very helpfulto my research. Her encouragement and support in various ways help mea lot. Her deep insights and passions in mathmatics inspire and enrich mygrowth as a student and a researcher.I gratefully thank the faculty and my fellow graduate students in theInstitute of Applied Mathematics. They provide a friendly, helpful, coor-perative environment for my research and study.I would also like to thank my parents. They have always been supportiveand encouraged me with their best wishes.xiDedicationTo my parentsxiiChapter 1Introduction1.1 BackgroundGenerally speaking, the term tipping point refers to an abrupt transition indynamical behavior observed as the system moves to a qualitatively differentstate, due to small changes in one or more factors [25]. Tipping pointshave been observed in many different fields, including the start and end ofice-ages [2], environmental regime shift [42, 43], synchronized behaviour ofneural activity [9], catastrophic collapse in ecology, and power systems [24].As many of these transitions are often irreversible or difficult to reverse,it becomes important to develop analytical techniques for predicting thelocation of tipping points and understanding the key contributing factors totipping.A number of analyses indicate that tipping points can depend on thebifurcation structure of the system, external noise or forcing, or the driftingrate of parameters that vary in time [4]. Of course in applications, a com-bination of these factors may contribute to the system simultaneously. Forexample, tipping can happen in the system with a saddle-node bifurcationpoint in which the stable and unstable equilibrium collide each other whenthe system transits from the stable state to unstable one and it’s irreversible.11.1. BackgroundA significant number of studies have focused on systems that have a slowlyvarying parameter approaching a saddle-node bifurcation point, where thetransition, also called tipping, occurs for parameter values beyond that ofthe static bifurcation point, known as a delayed bifurcation. Earlier recog-nition of this phenomenon is discussed in the context of optics [31] andneuronal dynamics [5], as well as general bifurcation theory [32]. Tippingcaused by the structure of a saddle-node bifurcation appears in many appli-cations. For example, the zero dimensional global energy balance model isdeveloped with a saddle-node bifurcation to explain the transition betweenthe stable global temperature and a ”deep freeze” situation [22], an energybalance model with the structure of a saddle-node bifurcation is explored toexplain the transition of the extent of Arctic sea ice between existence andextinguishment [3, 16], and a two-phase growth model with a saddle-nodebifurcaion is used for population resilience [24].An additional factor appearing in many applications is the periodic vari-ation. Often reduced models, such as those for the extent of Arctic sea ice[18], are developed for analysis by averaging over oscillations with higherfrequency, but a closer look at these types of models indicates different loca-tions of the tipping points when the oscillations are included. For example,the bifurcation point of a saddle-node bifurcation is shifted if the periodicoscillations are included. Seasonal forcing like the transimission rate orbirth/death rate changes the bifurcation diagram of the epidemic dynamicsystem[11], or induces the change between multistable behaviour and chaoticdynamics [30]. Periodic stimuli in nervous system triggers different bifur-cation structures on th mechanisms of the firing onset [36, 37]. Periodic21.1. Backgroundvariations can occur on different time scales, such as the seasonal variationconsidered in [16] or regular fluctuations on the time scale of decades [39].To explore the effect of periodic forcing with either low or high frequency,we consider the canonical model for a saddle-node bifurcation with a slowlyvarying bifurcation parameter with additive periodic forcing. In this settingwe explore how the location of the tipping point depends on the combinedeffect of the frequency, amplitude, and the drifting rate at which the bi-furcation parameter approaches the static bifurcation point. We use themethod of multiple scales and matched asymptotic expansions as in [6, 23]to analyze different relationships between the frequency of the oscillationand slow drifting rate, as well as the different ranges of the amplitude ofthe oscillation. In addition, to explore the influence of the amplitude on theposition of the tipping point, we introduce a concavity analysis to comparethe local behaviour of the system and its concavity threshold besides whichthe sign of the concavity is different, which essentially determines whethertipping happens.Many other studies focus on systems in which the drifting bifurcationparameter varies toward the saddle-node bifurcation point under the effectof an additive noise. The structure appears in many applications, includingmodels of global energy balance [22], extent of Arctic sea ice [3, 16], pop-ulation resilience [24], and regime shifts in gene expression dynamics [26],and so on. Focal points for prediction of tipping include statistical measuresbased on historical data that track the value of the slowly varying param-eter where an abrupt transition occurs [19]. These types of measures havebeen compared with canonical models, illustrating that in the presence of31.1. Backgroundnoise tipping can occur well before the slowly varying parameter reachesthe static bifurcation point of the model [38]. Whether the location of thetipping point is advanced, that is “early”, or delayed relative to the staticbifurcation point depends on the relationship between the drifting rate andthe strength of the noise, as studied for different underlying bifurcationstructures in [28, 34, 35]. Predictive techniques using statistical measuresand time series analyses have been developed to identify signatures of prox-imity to a tipping point or potential for early tipping, with examples in[19, 24, 41–43] and additional references in [27]. However, those measuresare not necessarily applicable if the tipping is driven by changes in the drift-ing rate or other factors, indicating that a better understanding of the roleof different tipping mechanisms is necessary.The concentration of the sample paths of the nonlinear system in theneighbourhood of the slowly varying equilibrium solution has been explored[28]. The delay of the transition caused by drifting bifurcation parameteris shown to be significantly reduced by small noise in many different bi-furcation structures [34]. The transition probability of a nonautonomousstochastic system with double well potentials away from the bifurcationpoint has been studied [35]. To detect the existence of a tipping point, dif-ferent indicators have been taken into account, such as critical slowing downof fluctuations indicating lost of resillience [24, 27, 41], and changes of theskewness [19, 42, 43]. Predictive techniques such as the method of movingwindows on time series can estimate the drifting rate and the amplitudeof the noise [38]. However, those techniques were built upon the canonicalmodels. The asymmetry of the bifurcation structure hasn’t been considered41.2. Literature Reviewand no explicit expression of the probability density function or probabil-ity of early tipping related to values of the bifurcation parameter has beengiven.To further understand the early tipping caused by the drifting bifurcationparameter and noise, we consider the canonical model for a saddle-nodebifurcation with a slowly varying bifurcation parameter and an additive noiseand explore the effect of the drifting rate and the amplitude of the noise byfinding the probability density function (PDF) of the system in a certainstate at a certain time. In order to capture the asymmetries of the PDF,instead of using Gaussian type approximation, we develop a new schemeof asymptotic analysis including patching different types of behaviours fordifferent parts of space and time domains. We use the WKB method and themethod of multiple scales to analyze different relationships of the driftingrate and the noise which can enhance the chance of an early tipping. To takeinto account the probability of the system staying in the neighbourhood ofthe unstable equilibrium, we consider a convection-diffusion process locallynear the unstable equilibrium to estimate the left tail of the PDF, which iscombined with the WKB approximation to provide the PDF of the systemand the escape probability. The approximation of the PDF can also be usedto analyze the indicators of early tipping.1.2 Literature ReviewFor the system with a saddle-node bifurcation structure and the bifurca-tion parameter slowly approaching the critical value, [32] finds the solutions51.2. Literature Reviewasymptotically in the neighbourhood of the stable and unstable branches.In the local region where the outer approximation is no longer valid, thesystem is in form of a Riccati equation and the solution is in form of theAiry function based on the exponential decreasing behaviour of the systemshortly after the bifurcation parameter crossing the critical value. Detailsare shown in Section 2.1.The effect of an external forcing on the saddle-node bifurcation with astatic bifurcation parameter has been studied [40]. Without a slowly varingbifurcation parameter, the fixed points on the stable and unstable bifur-cation branches are replaced by invariant curves. Similar results like theclassical results for the system with a saddle-node bifurcation are obtained.There exists a unique critical value for the bifurcation parameter. If thebifurcation parameter is smaller than the critical value, there exist two dis-tinct continuous invariant curves. If the bifurcation parameter equals tothe critical value, there exists only one invariant curve. If the bifurcationparameter is larger than the critical value, there is no invariant curve. Thebifurcation structure of the energy balance model for the extent of the Arc-tic sea ice described in [3, 15, 16] is similar to the one in [40]. If all theseasonal varying parameters are replaced by their annual average, the sys-tem has a double saddle-node bifurcation, with two bifurcation points, twostable branches connected by the one unstable branch. If the seasonal varia-tions such as the shortwave radiation flux are included, the annual change ofthe extent of the Arctic sea ice forms a bifurcation structure similar to theone in [40], given by numerical simulations based on constructing an annualPoincare´ map and finding a fixed point, with the bifurcation point slightly61.2. Literature Reviewshifted from the one based on the annual average. The bifurcation struc-ture of a Morris-Lecar neuron under a sinusoidal stimulus has been studied[36, 37]. A simple bifurcation structure with a bell-shaped firing onset isshown with the help of bifurcation and continuation algorithms. However,the time dependent bifurcation parameter which is necessary for the abrupt,irreversible transition is not considered. The influence of the frequency andthe amplitude of the periodic forcing hasn’t been explored.To count in all the factors, we first take the canonical model with asaddle-node bifurcation, including a slowly varying bifurcation parameterand an additive periodic variation. Due to the drifting rate of the bifur-cation parameter and the frequency of the oscillatory term, the method ofmultiple scales is required. Due to the relationship of the drifting rate ofthe bifurcation parameter and the frequency of the oscillation, we considerthe case with the high and low frequency separately. For the high frequencycase, we explore the relationship between the location of the tipping point,and the parameters from the periodic forcing. For the low frequency case,we introduce a concavity analysis in the local region compared with the localbehaviour of the system, to determine whether tipping happens at a certainvalue of the bifurcation parameter due to the amplitude threshold of theperiodic oscillation.For systems with a saddle-node bifurcation under the effect of an addi-tive noise, a lot of work has been done for the case with static bifurcationparameters. The Fokker-Planck equation (FPE), or the forward Kolmogorovequation is introduced to find the PDF of the system in a certain state ata given time [17]. With fixed bifurcation parameter, the WKB method is71.2. Literature Reviewused to find the so-called stationary distribution by solving the related FPEasymptotically [6, 17]. In case of drifting bifurcation parameters for a pitch-fork bifurcation and a Hopf bifurcation, for the noise much smaller than thedrifting rate of the bifurcation parameter, a Gaussian-type approximationof the PDF is given in [34] by local asymptotic analysis near the solutionof the deterministic system without the noise. However, for the noise com-parable to the drifting rate of the bifurcation parameter, the Gaussian-typeapproximation cannot capture the large bias on the left tail of the PDF,which indicates the possibility of the early tipping.To explore the probability of the early tipping before the bifurcationparameter reaches the critical value for the system with a saddle-node bifur-cation and a slowly varying bifurcation parameter, a sample-path approachwas introduced in [28]. The sample paths of the system are concentrated in aneighbourhood of the solution of the deterministic system. The relationshipbetween the amplitude of the noise and the drifting rate of the bifurcationparameter is given, showing that if the drifting rate of the bifurcation pa-rameter is slower than the square of the noise, the sample paths are likelyto escape from the neighbourhood of the stable branch before the bifurca-tion parameter reaches the critical value. However, no explicit expressionsfor the PDF are given. There are a lot of studies focusing on extractingthe nonlinear behaviour based on time series, assuming the time series isobtained from a system with a saddle-node bifurcation, a slowly varyingbifurcation parameter and an additive Gaussian noise [19, 38]. The loca-tion of the slowly drifting equilibrium is found by detrending, and the lineardecay rate, which is also used to find the amplitude of the noise, is found81.2. Literature Reviewby autocorrelation function [19]. The skewness and quasi-stationary densityare used as a sign of the nonlinearity based on the change of the potentialwell of the systems due to the decreasing distance between the bifurcationparameter and its critical value[19]. Even though in [42, 43], the trend ofthe skewness is claimed to be an indicator of the early tipping caused byregime shifts, the approximation of the skewness based on a time series in[19] shows that a significantly non-zero approximation of the skewness isobtained as a sign of nonlinearity but the change of the skewness based on atime series is negligible when the time series approaches tipping, comparedto the change of the skewness based on the empirical stationary density.However, in [19, 42, 43], the empirical stationary density of the system witha saddle-node bifurcation and an additive noise is computed based on a con-stant bifurcation parameter, which is related to a stationary process wherethe random realization is re-injected at +∞ whenever it has escaped to −∞.The effect of the slowly varying bifurcation parameter is not considered.To further understand the combined effect of the slowly varying bifur-cation parameter and the amplitude of the noise, an explicit expression ofthe probability density function for the system is necessary. It is time con-suming to get a precise approximation by the Monte-Carlo simulation andit is difficult to understand the parametric behaviour of the system behindthat. Hence, it is important to develop a general method to analyse theeffect of those factors and find their roles in the PDF, which in the end pro-vides the information of the early tipping. To capture the quasi-stationarybehaviour of the PDF relatively far from the bifurcation point, we use theWKB method on the related FPE and solve it asymptotically. To capture91.3. Summary of Resultsthe behaviour of the system near the unstable equilibrium of the saddle-node bifurcation when the bifurcation parameter is close to the bifurcationpoint, we combine the WKB method with the local analysis for a convection-diffusion process.1.3 Summary of ResultsIn Chapter 2, we use the method of multiple scales to find an explicit ex-pression of the location of the tipping point. For a high-frequency forcing,there are two main contributions to the location of the tipping point relativeto the static bifurcation point: a delay due to the slowly drifting bifurcationparameter, and an advance proportional to the square of the amplitude tofrequency ratio of the oscillation. The shift due to the high frequency pe-riodic forcing is also compared with the case of a white noise forcing. Forthe low-frequency forcing, we find that the location of the tipping point isaffected not only by the frequency, the amplitude, and the drifting rate, butalso by the behavior of the trajectory as the oscillations cause it to approachzero. This effect can be related to the phase of the oscillations, given by aninitial condition, and can lead to either delayed or advanced tipping. Wesee that unlike the high-frequency forcing, where the location of the tippingpoint varies smoothly with the other parameters, in the low frequency casethere are discontinuities in the location of the tipping point as a functionof, for example, the amplitude. For different relationships between the fre-quency and the drifting rate, we use a local analysis of the concavity of thetrajectory of the system to find critical values of the amplitude where the101.3. Summary of Resultstipping point location shifts dramatically.In Chapter 3, we develop a general method related to a WKB approx-imation and a convection-diffusion process to find the PDF of the systemwith a saddle-node bifurcation, a slowly varying bifurcation parameter andan additive white noise. To solve the FPE of the system in order to find thePDF, we introduce a WKB form of the solution as the value of the PDFrapidly decreases from the peak [6]. The WKB approximation is valid forthe bifurcation parameter relatively far from its critical value when thereis no existence of the early tipping. The leading order approximation isin form of the expression of the potential well of the system which is alsoknown as a quasi-stationary approximation. The approximation includinghigher order terms shows that the PDF either approaches to 0 or becomesa singularity at the unstable equilibrium which indicates that the WKB ap-proximation is no longer valid. For the bifurcation parameter closer to itscritical value where nonzero values of PDF rise near the unstable equilib-rium, the WKB approximation is not applicable in the neighbourhood ofthe unstable equilibrium, which requires us to consider the FPE separatelynear the stable and unstable equilibrium. We introduce a local expansionfor smaller bifurcation parameter and combine the PDF given by the WKBapproximation centered at the stable equilibrium and the PDF given by theconvection-diffusion process near the unstable equilibrium by the normaliza-tion of the whole PDF. With the approximation of the PDF at any value ofthe bifurcation parameter, we are able to explore the relationship betweenthe early tipping and the amplitude of the noise, the drifting rate of thebifurcation parameter, even the shape of the saddle-node bifurcation.111.3. Summary of ResultsIn Chapter 4, we use the general method developed in Chapter 2 andChapter 3 to analyse the location of the tipping point in a complex systemwith a saddle-node bifurcation, a slowly varying bifurcation parameter, anda periodic forcing or an additive noise. The Morris-Lecar model is consideredwith both high and low frequency oscillations separately to show that themethod of multiple scales can be applied in a system with higher dimensionsand the effects of the parameters from the oscillatory term on the location ofthe tipping points is similar to those in Chapter 2. The energy balance sea icemodel in [16] is analysed with a modified method based on the one in Chapter2 because the normalized amplitude of the periodic forcing is larger than thefrequency of the oscillation, which is different from the assumption for theanalysis in Chapter 2 for the high-frequency case. We apply the analysissimilar to those in Chapter 3 on the zero-dimensional climate model [22]with a saddle-node bifurcation, a slowly varying bifurcation parameter andan additive noise, to explore the effect of the slowly varying drifting rate,the amplitude of the noise, and the structure of the bifurcation diagramon the PDF of the system, and the probability of early tipping. For thebifurcation parameter relatively far from or closer to the critical value, ourmethod provides accurate approximation of the PDF, which also providesan approximation of the probability of the early tipping. The probabilityof the tipping shows that the asymmetry of the bifurcation diagram playsan important role to determine the probability of the early tipping, whichcannot be detected by the tipping indicators as skewness or variance of thePDF.In Chapter 2 and Section 4.2 we compare our analytical results to simu-121.3. Summary of Resultslations of the model using the second order Runge-Kutta method, while forthe Morris-Lecar model in Section 4.1, we use the fourth-order Runge-Kuttamethod. In Chapter 3 and Section 4.3, we compare our analytical results tosimulations of the stochastic model using the Euler-Maruyama method andthe Monte-Carlo simulation.13Chapter 2Tipping Points near aDelayed Saddle-nodeBifurcation with PeriodicForcing2.1 The ModelWe consider a canonical model with additive periodic forcing,dxdt = a− x2 +Asin(Ωt), dadt = −µ,x(0) = x0, a(0) = a0 > 0, (2.1)where a(t) denotes the slowly varying bifurcation parameter with the driftingrate 0 < µ ≪ 1, A is the amplitude of the oscillation (A > 0), and Ωrepresents the frequency. We assume a0, A are positive constants and x0 =√a0.The system (2.1) is generalized from the canonical model for a saddle-142.1. The Modelnode bifurcation when a(t) is a constant and A = 0 [13]. For any constanta > 0, there exist a stable equilibrium x+e =√a and an unstable equilibriumx−e = −√a. For a < 0, no equilibrium exists and the trajectory of (2.1)decreases exponentially for any initial condition x0. The bifurcation diagramis shown in Figure 2.1.For the system (2.1) without the periodic forcing term (A = 0), thereis no equilibrium because a(t) is slowly varying. With initial conditions x0and a0 > 0, the trajectory of the system is exponentially attracted to onespecific solution which we call the slowly varying equilibrium solution, ana-lyzed in [32] for A = 0 which we briefly summarize here. For a(t) > 0, theslowly varying equilibrium solution is close to x = √a for a constant whichwe call the stable branch. As a(t) crosses zero, there is a transition fromthe slowly varying equilibrium solution to the rapidly decreasing behaviour.This transition occurs for the value a(t) below the static saddle-node bifur-cation point a = 0, and thus is called a delayed bifurcation (see [32] andreferences therein). When such an abrupt qualitative change occurs, it is of-ten referred to as a tipping point regardless of the value of a(t). The methodof matched asymptotic expansions is applied to find the approximation ofthe slowly varying equilibrium solution for A = 0 [28, 32], using a slow timescale µt and solving for x(t) as a function of a(µt),x(a(µt)) =√a+ µ4a −532µ2a5/2 +O(µ3), (2.2)for a(µt) ∼ O(1). Notice that if a(µt) is O(µ2/3), (2.2) is no longer valid.152.2. High Frequency Oscillation Ω ≫ 1Then a local approximation is neededx(a(µt)) ∼ −µ1/3Ai′(a/µ2/3)Ai(a/µ2/3) , (2.3)where Ai denotes the Airy function. The asymptotic approximation of theslowly varying equilibrium solution of (2.1) with A = 0 is shown in Figure2.1, and compared with the static saddle-node bifurcation diagram. Thetipping happens near the singularity in the expression (2.3), which is thefirst zero of the Airy function,ad = µ2/3 · (−2.33810...) (2.4)The smaller that µ is, the slower a(t) varies, and the closer that the tippingpoint is to the static bifurcation point a = 0. The µ2/3 law characterizes thedelay of the jump transition and has been verified experimentally in variousphysical systems [12, 14, 21]. This scaling law reappears throughout theresults in this chapter.For A > 0 in (2.1), we explore the effect of the amplitude and frequencyof the oscillation on the location of the tipping point, considering both highand low frequency in combination with the slow drift of the parameter a(t).2.2 High Frequency Oscillation Ω ≫ 1In order to see how the tipping location is influenced by the high frequencyperiodic forcing, we first consider (2.1) with constant a (µ = 0) and de-termine the critical value of the bifurcation parameter where the bounded162.2. High Frequency Oscillation Ω ≫ 1−0.1 0 0.1 0.2 0.3−0.8−0.6−0.4−0.200.20.40.60.8axFigure 2.1: The asymptotic approximation of the slowly varying equilibriumsolution (2.2), (2.3) of the system (2.1) without the periodic forcing (A = 0)(black dash-dotted) is compared with the saddle-node bifurcation diagram(red solid and dashed lines). The solid (dashed) line is the stable (unstable)steady state. The solid circle indicates the static bifurcation point. Thevertical blue dotted line indicates the value of ad (2.4), the singularity of(2.3). The drifting rate is µ = 0.0025.attracting solution loses stability. We identify the bifurcation point denotedap, analogous to the saddle-node bifurcation point in Section 2.1, using themethods of multiple scales and matched asymptotic expansions. This resultshows that the periodic forcing shifts the location of the bifurcation pointto a value of a > 0. We compare this effect with the shift in the locationof the tipping point for combined slowly varying a(t) and high frequencyforcing. We obtain an expression for this shift in terms of the drifting rateµ of a(t), the frequency Ω, and the amplitude A. First we consider values of172.2. High Frequency Oscillation Ω ≫ 1A smaller than Ω or of the same order of magnitude of Ω. For larger valuesof A, a different approach is needed, as shown in Section 2.2.3.2.2.1 Constant bifurcation parameterFor µ = 0 (a(t) a constant) and A 6= 0, the solution of (2.1) oscillatesperiodically around√a for a an O(1) constant but it exponentially decreasesfor parameter values below the bifurcation point a = ap. Figure 2.2(LEFT)shows the bifurcation diagrams of (2.1) with µ = 0 and different periodicforcings. For larger amplitude A, the bifurcation point ap gets further fromthe critical value a = 0 for A = 0. The trajectories in Figure 2.2(RIGHT)show the stable and unstable periodic solutions overlap for larger A.We locate ap, using the method of multiple scales [6], based on both afast time scale T = Ωt and the original O(1) time scale t, the slow time scalein this case. For x = x(T, t), (2.1) becomesxT +Ω−1xt = Ω−1a− Ω−1x2 +Ω−1Asin(T ), (2.5)and the expansionx ∼ x0(T, t) + Ω−1x1(T, t) + Ω−2x2(T, t) + · · · , (2.6)is substituted into (2.5), yieldingO(1) : x0T = 0, (2.7)O(Ω−j) : xjT = Rj(T, t), j = 1, 2, 3 · · · . (2.8)182.2. High Frequency Oscillation Ω ≫ 10 0.05 0.1 0.15 0.2 0.25 0.3−0.4−0.200.20.40.6ax0 2 4 6 8 10−0.8−0.6−0.4−0.200.20.40.60.8txFigure 2.2: LEFT: Bifurcation diagrams for the average of x over one periodfor (2.1) with no drifting rate (µ = 0) and different periodic forcings: A = 0(Red); A = 3 (Black) and A = 5 (Blue). Frequency Ω = 4π. The solid(dashed) lines are the stable (unstable) steady states. The periodic solu-tions are generated by Poincare´ map as mentioned in [15]. By numericallyintegrating (2.1) based on the 2nd order Runge-Kutta method over one pe-riod, we can find the initial conditions that generate periodic solutions. Theblack (blue) triangle marks the numerical approximation of ap below whichthere exists no periodic solution. RIGHT: Numerical simulations of periodicsolutions of (2.1) with no drifting (µ = 0) and different periodic forcings.The black (blue) curve corresponds to the black squares (blue circles) inRIGHT for a = 0.095 and A = 3 (a = 0.15 and A = 5). The solid (dashed)lines are the stable (unstable) solutions. The trajectories are generated bythe 2nd order Runge-Kutta method starting from the initial conditions wefound in LEFT figure.192.2. High Frequency Oscillation Ω ≫ 1The leading order equation (2.7) indicates that x0 = x0(t) so that x0 dependsonly on the slow time variable t.To find xj for j = 0, 1, 2, · · · , we use a solvability condition to eliminatesecular terms in the solution of (2.8) that monotonously increase in the fasttime variable T . The solvability condition is that the right-hand side of (2.8)is orthogonal to the homogeneous solution of xjT = 0 [6],limT→∞1T∫ T0Rj(u, t)du = 0, j = 1, 2, 3 · · · . (2.9)Under the assumptions of the method of multiple scales, t and T are treatedas independent variables. For j = 1,R1(T, t) = a− x20 − x0t +Asin(T ). (2.10)By applying the solvability condition (2.9) on (2.10), we havex0t = a− x20, ⇒ x0 =√a,x1T = Asin(T ), ⇒ x1 = −Acos(T ) + v1(t),(2.11)where v1(t) is a function that needs to be determined in the higher orderanalysis.The higher order corrections xj , j = 2, 3 are obtained in Appendix A,yielding the asymptotic approximation for the attracting solution of (2.1)for µ = 0, a a constant, asx ∼√a+Ω−1 · [−Acos(T )] + Ω−2 ·[2√aAsin(T )− A24√a]+ · · · . (2.12)202.2. High Frequency Oscillation Ω ≫ 1Similarly, one can find the approximation of the solution for the case whenx0 = −√a and show it to be unstable.For A = 0, (2.12) reduces to x = √a with a saddle-node bifurcationpoint at a = 0. For A > 0, (2.12) is valid when both A and a are O(1). Fora = O(Ω−2), we substitute a = Ω−2b and introduce the expansion (2.6) into(2.5) to find ap. In Appendix A, we find the leading order approximation ofthe attracting solution for a = O(Ω−2)x ∼√a− A22Ω2 −AΩcos(T ) + · · · . (2.13)For a a constant, (2.13) shows that there are no attracting solutions if a <ap = A2/(2Ω2).The approximation of the critical value ap is no longer valid if A > Ωor Ω ≪ 1. Figure 2.3 shows the stable and unstable periodic solutions andtheir phase space at a > ap and a < ap where the frequency Ω is O(1) andthe amplitude of the oscillation A > Ω. We notice that for a < ap, shown asthe blue curves in Figure 2.3, the stable and unstable periodic solutions existand their trajectories overlap in phase space. In Section 2.2.3 and Section2.3, we explore the location of the tipping point for those cases under theinfluence of the slowly varying bifurcation parameter.2.2.2 High frequency oscillation in a slowly varying systemFor a(t) as in (2.1) with µ ≪ 1, we introduce a new parameter λ > 0by writing Ω = µ−λ in (2.1). The parameter λ captures the relationshipbetween frequency and drifting rate, which we relate to the tipping point212.2. High Frequency Oscillation Ω ≫ 10 2 4 6 8 10−505tx−5 0 5−6−4−20246810xxtFigure 2.3: LEFT: Numerical simulations of periodic solutions of (2.1) withno drifting (µ = 0) and the same periodic forcings with Ω = 2π and A = 10,but at different values of bifurcation parameter. Black: a = 8 > ap; Blue:a = 5 < ap. The solid (dashed) lines are the stable (unstable) solutions. Thetrajectories are generated by the 2nd order Runge-Kutta method startingfrom the initial conditions we found by Poincare´ map as mentioned in [15].RIGHT: Phase space for the stable and unstable periodic solutions in LEFT.222.2. High Frequency Oscillation Ω ≫ 1below. Here we consider the case A = o(Ω) for Ω ≫ 1, and we considerlarger values of A in Section 2.2.3. We again use the method of multiplescales, and find it is sufficient to use two scales, the fast time scale T = µ−λtand the slow time scale τ = µt. In terms of these time scales, (2.1) becomesµλ+1xτ + xT = µλ[a− x2 +Asin(T )], (2.14)aτ = −1.Then we focus on only (2.14) in the following analysis, analogous to [32]without periodic forcing. By substituting the general expansionx ∼ x0(T, τ) + µλx1(T, τ) + µ2λx2(T, τ) + · · · , (2.15)into (2.14) and solving for each order by applying the solvability condition(2.9) (details shown in Appendix B), we get the asymptotic approximationof x for a(τ) = O(1)x ∼√a+ µ4a + µλ [−Acos(T )] + · · · . (2.16)This result is similar to (2.2), with an oscillatory correction term µλ [−Acos(T )]from the high-frequency oscillatory term in (2.1). The asymptotic approx-imation (2.16) describes the attracting solution of (2.1) away from a = 0,but does not describe tipping. As a(τ) approaches zero, specifically fora(τ) ∼ O(µ2/3), (2.16) is no longer valid. In order to determine the influ-ence of the oscillations on the location of the tipping point compared to(2.4), a local analysis for a(τ) ∼ O(µ2/3) is needed.For the local analysis, we introduce a(t) = µ2/3α(s), the fast time vari-232.2. High Frequency Oscillation Ω ≫ 1able T = µ−λt and a new slow time variable s = µ1/3t and substitute into(2.1). The system (2.1) becomesµ1/3+λxs + xT = µλ+2/3α− µλx2 + µλAsin(T ), (2.17)αs = −1,and we focus on the local equation (2.17) to determine the location of thetipping point, analogous to the analysis in [32]. Separate analyses for dif-ferent values of λ suggest that it is useful to substitute x = −µλAcos(T ) +µ1/3y(T, s) in (2.17), yielding the asymptotic result for 13 ≤ λ < 1. Then wegetµλ+1/3ys + yT =µλ[µ1/3α− µ2λ−1/3A22 − µ2λ−1/3A22 cos(2T )+ 2µλAycos(T )− µ1/3y2], (2.18)and substitute the expansiony ∼ y0(T, s) + µλ+1/3y1(T, s) + · · · (2.19)into (2.18), yielding y0 = y0(s), and242.2. High Frequency Oscillation Ω ≫ 1O(µλ+1/3) : y1T =− y0s + α− µ2λ−2/3 ·A22 − µ2λ−2/3 · A22 cos(2T )+ µλ−1/3 · 2Aycos(T )− y20. (2.20)We apply the solvability condition (2.9) to (2.20), yieldingy0s = α− y20 − µ2λ−2/3A22 , ⇒ y0 = −Ai′(α− µ2λ−2/3 · A22 )Ai(α− µ2λ−2/3 · A22 ). (2.21)The local approximation of the attracting solution for a = O(µ2/3) is givenbyx ∼ µλ · [−Acos(Ωt)] + µ1/3 ·[−Ai′[(a− µ2λ A22 )/µ2/3]Ai[(a− µ2λA22 )/µ2/3]]+ · · · . (2.22)For 16 < λ < 13 , we substitute the expansion (2.15) into (2.17) and getthe same approximation as (2.22) in the end. For λ ≤ 16 , the frequency istreated as O(1) and is not considered here.Comparing (2.22) with (2.3), and recalling the relationship Ω = µ−λ, wefind that the position of the tipping point is determined by a that satisfiesAi[(a− A22Ω2 )/µ2/3] = 0. Then the tipping point with high frequency forcingand slowly drifting bifurcation parameter is given byahf ∼ ad + ap = µ2/3K +A22Ω2 , (2.23)252.2. High Frequency Oscillation Ω ≫ 1for K given in (2.4) and ap giving the same shift due to the periodic forcingas observed for the bifurcation point in (2.13) for µ = 0. Recall that in thissection we consider the case A = o(Ω), with larger values of A discussed inSection 2.2.3.−0.2 0 0.2 0.4 0.6−2−1.5−1−0.500.511.5ax−0.2 0 0.2 0.4 0.6−2−1.5−1−0.500.511.5axFigure 2.4: Tipping is shown for µ = .01. LEFT: The asymptotic approx-imations of the attracting solution for A = 1 and λ = .8, (2.16) (includinghigher order corrections up to O(µ3λ)) and (2.22) (blue dashed line) are con-sistent with the trajectory generated numerically (green dash-dotted line).The position of the tipping point is compared with the case of A = 0 (blackheavy dots). RIGHT: Blue dashed and green dash-dotted as on the LEFTwith A = 1 and λ = .2, compared with the numerical simulation for thestochastic process obtained by replacing the periodic forcing by an additivewhite noise (red solid line). ǫ = 0.2 is the amplitude of the noise for therealization shown. The numerical simulations are generated by the 2nd or-der Runge-Kutta method for the system (2.1) and by the Euler-Maruyamamethod for the stochastic process.Figure 2.4 shows the effect of the high-frequency oscillation in (2.1) onthe location of the tipping point. For higher frequency the location is almostthe same as the case with no oscillation but for lower frequency (smallervalues of λ) there exists the earlier tipping. We also see that a sloweroscillatory forcing in (2.1) results in an oscillation with larger amplitude262.2. High Frequency Oscillation Ω ≫ 1around the slowly varying equilibrium solution found without oscillatoryforcing (2.2), consistent with the expression for the outer approximation(2.16).In [28, 38], it was shown that the earlier tipping can be triggered alsoby white noise, in the system obtained by replacing the periodic forcing in(2.1) with an additive white noise, specifically, dx = (a − x2) dt + ǫdW (t)for W (t) a standard Brownian motion, which is also the canonical model(3.1) in Chapter 3. An earlier tipping is shown to occur with significantprobability for noise amplitude ǫ = O(µ1/2) or larger. Detailed analysis isshown in Chapter 3. For the shifts in the location of the tipping pointsshown in Figure 2.4 (RIGHT) resulting from noisy and periodic forcings,the amplitude of the oscillatory forcing (A = O(1)) is an order of magnitudelarger than the noise coefficient (ǫ = o(1)) that generates a comparableadvance of tipping. In the case of the high frequency periodic forcing, it isthe ratio of the amplitude to frequency that plays a significant role in thelocation of the tipping point, rather than the magnitude of the amplitudealone.The asymptotic approximation for the tipping point location ahf givenin (2.23) is compared with numerical simulations in Figure 2.5 for differentdrifting rates µ and varying exponents λ for Ω = µ−λ. There we see thecompetition between the two components in ahf . For A = 1 and smallervalues of µ, that is slower drifting rates, the location of the tipping pointis dominated by ad corresponding to A = 0 (2.4). In those cases only forsmaller values of λ, that is lower frequencies, is there a noticeable influenceof the periodic term, as seen in Figure 2.5(RIGHT). The advance of ahf due272.2. High Frequency Oscillation Ω ≫ 1to the periodic forcing as captured by ap clearly increases with A. For ex-ample, the case A = 5 is shown in Figure 2.5(LEFT), where there are largeradvances associated with ap for decreasing λ as compared with negligibleadvances for A = 1 and smaller λ. As noted above, for smaller values of λ(near λ ≈ 1/3 or below), the frequency Ω of the oscillatory forcing is notnecessarily large, depending on the value of µ. Then the expansion (2.15)is no longer a reasonable asymptotic approximation, and instead the behav-ior has some of the features discussed in Section 2.3 for the low frequencyforcing. As seen in Figure 2.5, the asymptotic approximation is valid for alarger range of λ for smaller values of µ. The approximation ahf and (2.22)are no longer valid for A ∼ O(Ω), and we consider results for larger valuesof A separately in the next section.−0.4 −0.2 0 0.2 0.400.20.40.60.81λa −0.4 −0.2 0 0.2 0.400.20.40.60.81λaFigure 2.5: The position of the tipping point for different frequencies Ω =µ−λ, shown for three drifting rates µ. At λ = 1 the values of µ showndecrease from left to right in both panels: µ = 0.05 in red; µ = 0.01 inblue; µ = 0.001 in black. The asymptotic approximation (2.23) is shown bythe solid line. The dash-dotted line shows the location of the tipping pointobtained from numerical simulations, defined as the value of a for x = −10.LEFT: A = 5; the dotted vertical lines show the position of the tipping pointad for A = 0; RIGHT: A = 1; the black dotted line shows the approximationof the bifurcation point ap = A2/(2Ω2) for µ = 0.282.2. High Frequency Oscillation Ω ≫ 12.2.3 High frequency oscillation with large amplitude(A≫ 1)In Section 2.2.1 and 2.2.2, the asymptotic approximations (2.16) and (2.22)are not valid if A = O(Ω) or larger. However, for some larger values A andhigh frequency Ω, we can rescale (2.1) to obtain a similar system forced byan oscillation with amplitude A = 1. In certain cases we can apply theapproximations from either Section 2.2.2 or Section 2.3.1 to the rescaledsystem.Substituting x =√Az and t = S/√A into (2.1), we getdzdS = h− z2 + sin(ωS),dhdS = −M,v(0) = v0 =x0√A, h(0) = h0 =a0A , (2.24)where ω = Ω/√A, h = aA , and M = A−3/2µ are the scaled frequency,bifurcation parameter, and drifting rate, respectively. In order to considerthe case A ≫ 1, we assume A = ΩP with P ≥ 1. Then the relationshipbetween ω and M is given by ω = Mγ , with γ a function of λ and Pγ = P − 23P + 2/λ. (2.25)For 1 ≤ P < 43 , the high frequency approximation (2.22) can be appliedto (2.24) if γ < −16 , which indicates that λ in Ω = µ−λ satisfies λ >212−9P ≥ 23 . For larger values of A (P ≥ 43), the oscillatory forcing and292.2. High Frequency Oscillation Ω ≫ 1the drifting bifurcation parameter both dominate the dynamics to leadingorder, so that the asymptotic approximation (2.16) dominated by the slowlydrifting parameter a(t) is no longer valid.For P > 2 and λ > 0, γ is positive which indicates that even though (2.1)has a high-frequency oscillatory forcing, in the normalized system (2.24) therelationship between ω andM corresponds to a low-frequency case, discussedin the next section. We note that γ is a monotone increasing function withrespect to P with a horizontal asymptote at γ = 13 .Summary of results: For a high-frequency forcing, the trajectory is ap-proximated by the linear combination of the slowly varying equilibrium so-lution found for A = 0 in [32] and an additive contribution that takes theform of the oscillatory forcing with a scale factor. Near the tipping point,the trajectory can be expressed in terms of an Airy function on the slowtime scale, with the oscillations averaged on this time scale shifting the sin-gularity that corresponds to the tipping point. The resulting expressions forthese contributions then give an expression for the location of the tippingpoint that varies continuously with the parameters. There are two compet-ing contributions to the location of the tipping point relative to the staticbifurcation point at a = 0. The first contribution is a delay due to theslow drifting rate µ of the bifurcation parameter a, well-known in a varietyof studies without periodic forcing [31, 32]. The second contribution is anadvance proportional to the square of the amplitude to frequency ratio ofthe oscillation, observed also in the case where the bifurcation parameter isconstant (µ = 0). Whether the combined effect results in an advance or de-lay compared to a = 0 depends on the magnitude of µ2/3 relative to A2/Ω2,302.3. Low Frequency Oscillation Ω ≪ 1showing the influence of not only the drifting rate µ of the bifurcation pa-rameter, but also the oscillatory forcing. In particular, we see noticeableadvances of the tipping point for combinations of larger amplitudes withfrequencies that are not very large. However, for the case A = O(Ω) orlarger, it is necessary to rescale the system (2.1) to get a similar one with anoscillation with amplitude A = 1. By comparing the normalized driftng rateM and frequency ω, the system can be either in the high or low frequencycase.2.3 Low Frequency Oscillation Ω ≪ 1For the case of a low-frequency oscillation (Ω ≪ 1), we introduce a newparameter ν > 0 by writing Ω = µν , capturing the relationship between thefrequency and the drifting rate. In the following we consider two differentcases: ν ≥ 1 so that Ω = O(µ) or smaller, or Ω = µν with 0 < ν < 1.2.3.1 Ω = O(µ) or smallerHere it is convenient to write Ω = c·µ for c a positive constant and to rescaletime as τ = µt as in [32]. Then we consider x as a function of a(τ) = a0− τin (2.1), to get− µdxda = f(a)− x2 = a− x2 +Asin(c · (a0 − a)). (2.26)312.3. Low Frequency Oscillation Ω ≪ 1Substituting x ∼ x0 + µx1 + µ2x2 + · · · into (2.26), we find the first twoterms of the approximate solution as a function of a(τ)x ∼ ±√f(a) + µ · f′(a)4f(a) + · · ·∼ ±√a+Asin(c · (a0 − a)) + µ ·1− c · Acos(c · (a0 − a))4(a+Asin(c · (a0 − a)))+ · · · , (2.27)where the positive sign corresponds to the solution that is attracting. Notethat the oscillatory term appears in the leading order contribution in (2.27),in contrast to the leading order approximation in the case of the high fre-quency oscillations (2.16) where there is no oscillatory term. Similar to (2.2),(2.27) is no longer valid for f(a) ∼ O(µ2/3). As in the cases considered inprevious sections, we expect a tipping point near a = ar, where f(a) > 0for a > ar and f(ar) = 0. A local analysis near a = ar is used to determinethe location of the tipping point. In our analysis and simulations we takea0 > ar for f(ar) = 0 with |a0 − ar| = O(1), ensuring that the system (2.1)attracts to the outer solution (2.27) prior to any tipping.Figure 2.6 shows numerical simulations of (2.1) for different values ofA and fixed a0, illustrating important characteristics of tipping that arethe basis for our analysis. The heavy black dots correspond to A = 0 (nomodulations). For larger values of A, the trajectory x = x(a) exhibits largeroscillations before tipping on its approach to the static bifurcation pointat a = 0. As could be expected, the value of a where the tipping occursincreases with A for a > 0. For amplitudes A where the trajectory crossesx = 0 for a(t) > 0, tipping can occur well before the static bifurcation pointis reached. Abrupt shifts in tipping point location are observed for critical322.3. Low Frequency Oscillation Ω ≪ 10 5 10 15 20−3−2−10123456ax (I)(II)(III)(IV)(V)−2 0 2 4 6 8 1012345678910AaFigure 2.6: Parameter values are µ = Ω = 0.01(c = 1), a0 = 20 for allgraphs. LEFT: Numerical simulation of (2.1) by the 2nd order Runge-Kuttamethod. The Roman numerals correspond to cases with different values ofA: Case (I), A = 8.9647; Case (II), A = 8.9797; Case (III), A = 2.5268;Case (IV), A = 2.5418; Case (V), A = 5. With the exception of Case (V),A is near the critical value Ac at which there is an abrupt change of thetipping position. Heavy black dots: A = 0. The solid black square (redcircle) corresponds to parameter values (a,A) = (am, Am), in Cases (I) and(II) ((III) and (IV)). Note that Am is slightly less than the value of A forCases (I)-(IV). The dashed rectangle indicates the region shown on a largerscale in Figure 2.7(LEFT). The blue triangle marks the value of ar wheref(ar) = 0 for A = 5. RIGHT: The location of the tipping point as a functionof A in (2.1) obtained by numerical simulations, with tipping defined as thevalue of a at x = −10 as in the previous cases. Different markers indicatethe value of A and corresponding tipping point value, for the trajectoriesshown in the LEFT figure. Solid (black) square: (I); Open (blue) square:(II); Solid (red) circle: (III); Open (red) circle: (IV); Solid (blue) triangle:(V).332.3. Low Frequency Oscillation Ω ≪ 1amplitudes denoted Ac that are slightly greater than Am, where the pair(a,A) = (am, Am) satisfies f(am) = f ′(am) = 0 as discussed below (2.31).For example, the Cases (I) and (II) shown in Figure 2.6(LEFT) have thevalues of A = 8.9647 and A = 8.9797, respectively, with both values slightlygreater than Am ≈ 8.9485 for a0 = 20. However, the tipping for Case(I) occurs at a ≈ 3.80, (solid black square in Figure 2.6(RIGHT)), whilethe tipping for Case (II) occurs at a ≈ 8.82 (open blue square in Figure2.6(RIGHT)). Here the value of am is indicated by the solid black square inFigure 2.6(LEFT), which is near the tipping point for Case (II) but not forCase (I).Figure 2.6(RIGHT) shows the location of the tipping points for all tra-jectories represented in Figure 2.6(LEFT), and indicates the critical valuesAc where there are abrupt changes in the tipping point as a function of A.As we study separately below, these values Ac are near Am, and are relatedto the local behavior of the trajectory near am. Figure 2.7(LEFT) zoomsin on the trajectories (I) and (II) near a = 8.9 and illustrates the dramaticchange of concavity near am.Before we investigate the solution near (am, Am) where f(am) = f ′(am) =0 (and near Ac), we first consider the case where tipping occurs close to a =ar where f(ar) = 0 and f ′(ar) = O(1), as in case (V) in Figure 2.6(LEFT).Near a = ar, specifically f(a) = O(µ2/3), we substitute a− ar = µ2/3B andx = µ1/3X into (2.26), yieldingO(µ2/3) : −XB = f ′(ar)B −X2 = [1− cAcos(c · (a0 − ar))]B −X2, (2.28)342.3. Low Frequency Oscillation Ω ≪ 1which has the solution for f ′(ar) 6= 0,X = −C1/3Ai′(C1/3B)Ai(C1/3B) ,⇒x = −(µC)1/3Ai′((C/µ2)1/3(a− ar))Ai((C/µ2)1/3(a− ar)), C = f ′(ar). (2.29)Similar to the analysis for (2.3) and (2.22), the location of the tippingpoint is determined by the singularity of (2.29) where Ai((C/µ2)1/3(a −ar)) = 0, which isalf = ar + ad/(f ′(ar))1/3 for f(ar) = 0 and f(a) > 0 for a > ar, (2.30)where ad is given in (2.4). We differentiate f(ar) = 0 with respect to A, toobtain dar/dA = (CA)−1ar, which indicates that the location of the tippingpoint increases (decreases) with increasing A for ar negative (positive), asshown in the Figures 2.6(RIGHT) and 2.9(RIGHT).Noting that (2.29) is obtained from (2.28) for f ′(ar) 6= 0, we propose aseparate analysis for (a,A) close to (am, Am), which satisfies the conditions,f(am) = am +Amsin(c · (a0 − am)) = 0, (2.31)f ′(am) = 1− cAmcos(c · (a0 − am)) = 0, ⇒ c2a2m + 1 = c2A2m. (2.32)The behavior of f(a) near am for A = Am is shown in Figure 2.7(LEFT).Figure 2.7(RIGHT) shows Am as a function of am for different values of a0as in (2.31) and (2.32). From these conditions, we see that for each a0 there352.3. Low Frequency Oscillation Ω ≪ 1is a family of initial values a0 + 2kπc for k an integer that corresponds to thesame pair (am, Am). Equations (2.31) and (2.32) indicate that Am > 1/cand 0 < am < Am. We take A > 0 in (2.1), noting that similar results canbe found for A < 0.8.7 8.8 8.9 9 9.1 9.2−1−0.500.51a (I)(II)0 2 4 6 8 1012345678910aA(I)&(II)(III)&(IV)Figure 2.7: For all graphs, µ = Ω = 0.01 (c = 1). LEFT: Figure 2.6(LEFT)zoomed-in near am for f ′(am) = 0 (black square), with f(a) for a0 = 20 andA = Am = 8.9485 given by the red dashed curve. Numerical simulations ofCase (I) (solid black line, A = 8.9647) and Case (II) (dash-dotted blue line,A = 8.9797) are also shown. RIGHT: Am as a function of am (solid red line),with a0 varying along this line as in (2.31) and (2.32). The black square(redcircle) marks the value of am and Am for Cases (I) and (II) (Cases (III) and(IV)), with a0 = 20 as in Figure 2.6.A comparison of cases (I) and (II) in Figure 2.6(LEFT) with Figure2.7(LEFT) suggests that the critical amplitude Ac, for which there is a jumpin the location of the tipping point, is near Am, with Ac > Am. Further-more, the trajectories shown in Figure 2.7(LEFT) have different concavitiesnear a = am for A above and below Ac. Therefore, to find Ac and the corre-sponding tipping point, we consider the concavity of x in (2.26) with respect362.3. Low Frequency Oscillation Ω ≪ 1to a near am. Using (2.26), we write x′′(a) in terms of the local variablesa− am = µη, x = µξ, and A = Am + µA1 , (2.33)to focus on the neighbourhood of a = am and on amplitude values near Am.Substituting the local variables (2.33) in (2.31) and (2.32), and expandingsin(cµη) and 1− cos(cµη) for µ≪ 1, we obtain the leading order expressionfor the local concavity of ξ(η),d2ξdη2 = µ( A1Am− c2amη +2A1amAmξ)+O(µ2). (2.34)From (2.34) we find that the leading order expression of ξ = x/µ thatseparates trajectories of different concavities is given byξc =c2Am2A1η − 12am. (2.35)Then for ξ > ξc (ξ < ξc), near a = am the trajectory is concave up (down).To determine the critical value A1 (and thus Ac) for which the trajectorieschange concavity we compare (2.35) to the trajectory itself in terms of thelocal variables η, ξ, and A1. The local expression for the trajectory ξ isfound by substituting (2.33) into (2.26), and using the expansion ξ ∼ ξ0 +µξ1 + µ2ξ2 + · · · ,ξ ∼ A1amAmη + c0 + µ · [−(12c2am −A21a2mA2m)η33+ A1Am(1 + 2c0am)η22 + c20η + c1] + · · · , (2.36)372.3. Low Frequency Oscillation Ω ≪ 1where c0 and c1 are undetermined constants. It is the comparison of theslope of the concavity threshold (2.35), c2Am/(2A1), to that of the trajectory(2.36), A1am/Am, that yields the result for A1. If A1am/Am < c2Am/(2A1)for a − am = O(µ) and x = O(µ), then trajectories that enter the concaveup region stay there, as shown in Figure 2.8(LEFT), and tipping does notoccur near η = 0 (a = am). Alternatively, if A1am/Am > c2Am/(2A1),then locally the trajectory has a slope greater than that of the concavitythreshold and trajectories that are near or below (2.35) for a near am crossinto or continue in the concave down region, resulting in tipping as shown inFigure 2.8(RIGHT). Then the critical value of A1 for the change in concavitysatisfiesA21 =Ωc2A2m2µ√c2A2m − 1, c = Ωµ ⇒ Ac ∼ Am + µA1, (2.37)where we have used (2.32) to eliminate am and (2.33) to get Ac.To arrive at (2.37), we assumed that for A < Ac, the trajectory (2.36) isabove the concavity threshold (2.35) in the local region (2.33). To verify thisassumption, we compare the ξ-intercept of the concavity threshold (2.35),− 12am , to that of the trajectory, which to leading order is c0. Then c0 isdetermined by matching the local expansion (2.36) with the outer expansion(2.27), as shown in Appendix C,c0 = −12am+ C1 + C2.The ξ-intercept of the trajectory, relative to −1/(2am) is obtained from the382.3. Low Frequency Oscillation Ω ≪ 1sign of C1+C2. As discussed in Appendix C, for A1amAm <c2Am2A1 (A1amAm >c2Am2A1) C1+C2 > 0 (C1+C2 < 0), which is consistent with our assumptions in thederivation of (2.37). Figure 2.8 illustrates the local behaviour by comparingthe trajectory and the concavity threshold near a = am.8.7 8.8 8.9 9 9.1−0.500.5a 8.7 8.8 8.9 9 9.1−0.500.5aFigure 2.8: View of Cases (I) and (II) as in Figure 2.6(LEFT), zoomed in onthe scale |a− am| = µ1/3, with |a− am| = µ indicated by red dashed lines.The black square is the value of am where f(am) = f ′(am) = 0. LEFT:The trajectory of case (I) (solid black line) is compared with the concavitythreshold (2.35) (black dashed line with markers). RIGHT: The trajectoryof case (II) (blue dash-dotted line) is compared with the concavity threshold(2.35) (blue dashed with markers).The panels in Figure 2.9 illustrate the dependence of these results on theratio c = Ω/µ, and compare the tipping point locations obtained by asymp-totic approximations (2.30) and (2.37) and by numerical simulations. Thevariation of Am and Ac with the parameter c is observed by comparing Fig-ure 2.9(LEFT) and Figure 2.7 (RIGHT), and noting the values of Ac wherethere is a jump in the location of the tipping point in Figure 2.9(RIGHT).There we see that these jumps occur near A = 9 and near A = 2 for thecases with c = 1 (blue and green lines), in contrast to the cases with c = .2392.3. Low Frequency Oscillation Ω ≪ 10 2 4 6 8 105678910aA−10 −5 0 5 1012345678910AaFigure 2.9: LEFT: Values of Am (red solid line) from (2.31) and (2.32) as afunction of am for c = 0.2. The black square marks the value of am and Amnear parameter values where there is an abrupt change in the tipping pointlocation, shown for the cases with c = .2 in RIGHT. RIGHT: The positionof the tipping point as a function of amplitude A, comparing the asymptoticapproximation (2.30) with jumps at values (2.37) (solid lines) with numericalsimulation, defining the tipping point as the value of a when x = −10 (dash-dotted lines with markers). Different colours indicate different combinationsof frequency Ω and drifting rate µ. Blue: Ω = 0.01, µ = 0.01, and a0 = 20;Red: Ω = 0.01, µ = 0.05, and a0 = 32; Green: Ω = 0.1, µ = 0.1, anda0 = 20; Black: Ω = 0.1 and µ = 0.5, and a0 = 32.402.3. Low Frequency Oscillation Ω ≪ 1where 6 < Ac < 7 (red and black lines). Figure 2.9 also shows good agree-ment between results for the location of the tipping point as a function ofA obtained from simulations and from the asymptotic approximation forΩ ≪ 1.2.3.2 Ω = µν (0 < ν < 1)For c much larger than 1, so that the low frequency is much larger thanthe drifting rate, µ ≪ Ω ≪ 1, the approximation (2.37) is not valid. Nev-ertheless, simulation results shown in Figure 2.10(LEFT) for Ω relativelysmall illustrate phenomena similar to that shown in Figure 2.9(RIGHT).For a given a0 there are values Ac at which there are abrupt changes in thelocation of the tipping point. As in the previous case of Ω = O(µ), tippingis related to a change of concavity in the trajectory near a local minimumof f(a) at a = amin. However, for larger values of Ω, f(amin) < 0 may beO(1) for amin near the tipping point, in contrast to f(am) in Section 2.3.1which vanishes as shown in Figure 2.7.To approximate Ac for Ω = µν , we first consider the outer approximationto (2.26) given by (2.27) with c = µν−1. It is again valid only for f(a) > 0,since as f(a) approaches 0, the terms shown in (2.27) are both of the sameorder for µ ≪ 1. To estimate Ac in this case, we must consider multiple412.3. Low Frequency Oscillation Ω ≪ 1local minima amin of f(a) for a given a0, satisfyingf(amin) = amin +Asin(µν−1(a0 − amin)) < 0, (2.38)f ′(amin) = 1− µν−1Acos(µν−1(a0 − amin)) = 0,⇒ cos(µν−1(a0 − amin)) = µ1−ν1A. (2.39)Since we are considering larger values of the frequency Ω, the values of aminare relatively close together and close to values where f(a) = 0. Then wecan not use a local analysis near values of a where f(a) = 0, as was donenear ar in (2.28), but instead restrict our analysis near amin for f ′(amin) = 0.Substituting a = amin + µη and x = µξ into x′′(a) using (2.26), (2.38),and (2.39), we get the local expressions for the trajectory and its concavity.Keeping the terms that dominate these expressions for µ≪ 1 yieldsµdξdη = −f(amin), (2.40)µd2ξdη2 = µ2ν(f(amin)− amin)η − 2µf(amin)ξ. (2.41)As in Section 2.3.1, we look for the critical value of A = Ac that correspondsto the value where the local slope of the trajectory is equal to that of theconcavity threshold. The latter is obtained from setting ξ′′(η) = 0 in (2.41)yieldingξc =µ2ν(f(amin)− amin)η2µf(amin). (2.42)Equating the coefficient of η from the right hand side of (2.42) with ξ′(η) in422.3. Low Frequency Oscillation Ω ≪ 1(2.40) gives the relation−2f2(amin) = µ2ν(f(amin)− amin). (2.43)Then, using (2.38)-(2.39) and (2.43), we can eliminate f(amin) and obtain asystem of equations for the pair (amin, Ac),amin =√A2c −µ2Ω2 −Ω√2[A2c −µ2Ω2]1/4, (2.44)amin = a0 −µΩ[cos−1( µAcΩ)− 2kπ], (2.45)where k is a positive integer. Note that for fixed a0 the system of equations(2.44)-(2.45) is valid for µ = o(Ω). It has more than one solution, that is,more than one pair (amin, Ac) at which there is a change in concavity andthus a jump in the location of the tipping point as a function of A. Then forsmall µ or sufficiently large Ac, cos−1(µ/(AcΩ)) ∼ π/2 in (2.45), so that thejumps in the value of amin are dominated by 2kπµ/Ω for integers k, relatedto the multi-valued function cos−1(·).Figure 2.10(LEFT) compares the position of the tipping points obtainedfrom numerical simulations with the asymptotic approximation for the crit-ical values amin and Ac given by (2.44)-(2.45). Figure 2.10 illustrates thisapproximation for µ = .1, both for Ω = µν for ν ≈ .3, and for Ω = O(1).While the latter case does not appear to be covered by the asymptotic anal-ysis, we recall that in Section 2.2.3, we indicated that the case A > Ω2 withΩ = µ−λ for λ > 0 can be studied as a low frequency case by rescaling (2.1)into the normalized system (2.24). This rescaled system, with frequency432.3. Low Frequency Oscillation Ω ≪ 12 3 4 5 6 7 8345678910Aa 0 5 10 15 20−4−20246axFigure 2.10: The drifting parameter is µ = .1 for all graphs.LEFT:Numerical simulations of the position of the tipping point (dash-dotted lines), defined as the value of a when x = −10, compared withmarkers indicating (amin, Ac) from (2.44) and (2.45), the critical values ofthe amplitude at which there is an abrupt change in the location of the tip-ping point. Red dash-dotted line and solid squares are for the case Ω = 0.5and a0 = 20. Black dash-dotted lines and solid circles are for the caseΩ = 1.2 and a0 = 20. RIGHT: Numerical simulations of (2.1) with A nearAc in the location of the tipping points in (LEFT): A = 8.7 (red solid) andA = 8.9 (red dash-dotted), both with Ω = 0.5; A = 7.1 (Black solid) andA = 7.3 (Blue dash-dotted), both with Ω = 1.2. The black square and redcircle mark the related amin for Ω = .5, 1.2, respectively, satisfying (2.44)and (2.45).ω < 1 and slow drift parameter M ≪ 1, can be analyzed using expressionsanalogous to (2.44)-(2.45). For µΩ ≪ 1, we see that (2.44)-(2.45) providea good approximation for the location of the tipping point and the valuesof A for which there is a jump in location, even though we do not have ananalytical approximation for the solution x for Ω = O(1) and µ ≪ 1. Forlarger values of Ω with µ fixed, there are more values of Ac where this jumpoccurs, and the corresponding values of a = amin are closer together. Theincreasing number and reduced distance between these jumps indicate thatas Ω increases further, the relationship between the amplitude A and the442.3. Low Frequency Oscillation Ω ≪ 1location of the tipping point is a continuous function, as shown in Section2.2.2 for the large frequency oscillation.−2 0 2 4 6 8 100246810AaFigure 2.11: Position of the tipping point as a function of A in (2.1) obtainedby numerical simulation, defined as the value of a where x = −10 ( dash-dotted line with solid circles) compared with the asymptotic approximation(solid line) for Ω = 0.1(red); Ω = 5(blue); and Ω = 10(magenta) withµ = 0.1 and a0 = 20. The black dash-dotted lines give the location of thetipping point determined numerically for Ω = 1.2, which are compared withthe black circles indicating the values of A where there are jumps in thelocation of the tipping point based on the asymptotic approximation (2.44)and (2.45).An analysis of large amplitude oscillatory forcing also suggests a rescalingof the system, which allows the asymptotic approximations to be extended tosome cases where the frequency is not particularly large or small relative tothe drift rate of the bifurcation parameter. In Figure 2.11 we compare resultsfor low, high, and O(1) frequencies Ω, illustrating how the the location of452.3. Low Frequency Oscillation Ω ≪ 1the tipping point varies with the amplitude and frequency of the periodicforcing.Summary of results: For Ω = O(µ) ≪ 1, the trajectory is dominated bya nonlinear combination of the slowly drifting bifurcation parameter and theoscillatory forcing, since they give contributions that are of the same orderof magnitude. The location of the tipping point is given by the singular-ity of an Airy function, depending not only the amplitude, frequency, anddrifting rate, but also on the phase of the oscillation. The analytical expres-sion indicates that tipping occurs near values of a where the trajectory isnear x ∼ 0. The values of the parameter such as the amplitude where theanalytical expression (2.30) is no longer valid, captures the discontinuous de-pendence of the tipping location. The change of concavity for x ∼ 0 providesan analytical expression for the amplitudes showing discontinuities. Thesediscontinuities follow directly from the fact that tipping occurs only for val-ues of a where the trajectory is near zero: a slight change of the amplitudecan change the approach to zero of the trajectory, postponing tipping to alater point at which the trajectory crosses zero. For µ≪ Ω ≪ 1, the criticalvalues of the amplitude for the discontinuities are obtained by a similar lo-cal analysis for the change of the concavity. The shift in the location of thetipping point is on the order of one period of the low frequency oscillation.The discontinuity of the location of the tipping point due to the parame-ters of the oscillation has not been analysed in previous studies [36, 37, 40],which focus on the systems with oscillatory forcing and constant bifurcationparameter .46Chapter 3Tipping Points near aDelayed Saddle-nodeBifurcation with AdditiveNoise3.1 Scaling Relationship of the Drifting Rate andthe Amplitude of the NoiseSimilar to (2.1), we consider a canonical model with an additive white noise,dx = (a− x2)dt+ ǫdWt,dadt = −µ, (3.1)x(0) = x0, a(0) = a0 > 0,where ǫ denotes the amplitude of the additive noise and Wt represents thestandard Wiener process. For ǫ = 0 and µ = 0, (3.1) has a saddle-nodebifurcation at (ae = 0, xe = 0), which is the same as the one shown inFigure 2.1.473.1. Scaling RelationshipLet p(x, t) be the PDF of x at time t. The Fokker-Planck Equation(FPE) corresponding to (3.1) ispt = −[(a− x2) · p]x +ǫ22 · pxx,p(x, 0) = p0(x).(3.2)To capture the influence of the relationship of the slow drifting rate andthe amplitude of the noise on the process x, we introduce a new parameterβ > 0 by writing µ = ǫβ. In order to consider regions far from and close tothe bifurcation point (ae = 0, xe = 0), we introduce scaled variablesa = a0 − µt = ǫαs, x = ǫα/2ξ, (3.3)where αmeasures the distance between a and its critical value 0. x is rescaledto O(ǫα/2) in (3.3) because x ∼ √a0 − µt following the upper branch of thebifurcation diagram in Figure 2.1. Then (3.2) becomes− ps = ǫ32α−βp · 2ξ + ǫ 32α−βpξ(ξ2 − s) + ǫ2−βpξξ2 . (3.4)We assume that p(ξ, s) has the WKB form of the solution [6]p(ξ, s) ∼ exp{ψ(ξ, s)/ǫγ}, (3.5)where γ is non-negative and determined by α and β. Substituting (3.5) into(3.4), we get483.1. Scaling Relationship− ψs = ǫ32α−β+γ · 2ξ + ǫ 32α−βψξ(ξ2 − s) + ǫ2−β−γψ2ξ2 + ǫ2−βψξξ2 . (3.6)To find a non-trivial solution, we need at least two terms in the leadingorder. By comparing the O(ǫ 32α−β), O(ǫ2−β−γ), and O(1) terms, we findthree cases based on different values of α and β.For α < 23β, α < 43 and γ = 2 − 32α, the leading order O(ǫ32α−β) termsin (3.6) satisfy0 = −(s− ξ2)ψξ +12ψ2ξ . (3.7)Solving (3.7) yields a quasi-stationary approximation and the PDF isp(x, t) ∼ exp{ 1ǫ2 [−2(−ax+x33 ) + F0(a)]}. (3.8)For µ = 0, F0(a) is a constant and (3.8) is the stationary solution of theFPE for (3.1) with F0(a) = 0 [17]. For µ > 0, F0(a) and a are slowly varyingfunctions of the time scale τ , with F0(a) determined together with higherorder contributions as shown below.The approximation (3.8) reflects the effect of the noise from the aspectof the potential well. As a decreases, the depth of the left barrier of the welldecreases as shown in Figure 3.1. A small random forcing is likely to push(3.1) to the unstable state before the barrier disappears, which triggers theearly tipping. The relationship between the drifting rate of the bifurcationparameter and the amplitude of the noise given above is consistent with493.1. Scaling Relationshipthe results corresponding to a sample paths method introduced in [28, 38].It has been shown that if the nosie ǫ is larger than the square root of thedrifting rate of the bifurcation parameter µ1/2, the sample paths are likelyto escape from the neighbourhood of the stable branch for positive a of orderǫ4/3.−1 −0.5 0 0.5 1−0.4−0.3−0.2−0.100.10.20.30.4xU a=0.2a=0.01Figure 3.1: Potential well of system (3.1) at a = 0.2 (blue dash-dotted) anda = 0.01 (red solid). Solid (Open) circles represent the stable (unstable)equilibrium of (3.1) with µ = ǫ = 0.For α > 23β, β < 2, and γ = 2− β, the leading order of (3.6) are− ψs =12ψ2ξ . (3.9)Then a diffusion equation governs the leading order behavior, and a Gaussian-type approximation for the PDF is appropriate as in [34],503.2. WKB Approximationp(x, t) ∼ C√g(t)ǫ exp{−[x− F (t)− 2ǫ√gn]2· g(t)2ǫ2}, (3.10)where C is the normalization constant, F (t) is the solution of the determin-istic system (ǫ = 0 in (3.1)), g(t) and gn(t) are functions related to F (t)(details are shown in Appendix E). The approximation (3.10) is similar tothe Gaussian distribution, except the mean and variance are t-dependent,which indicates that (3.1) is centered on the slowly varying deterministicdynamics described by F (t). This relationship is also consistent with theresults obtained based on the sample paths method [28], which says that ifthe amplitude of the noise ǫ is smaller than the square root of the driftingrate µ1/2, the sample paths are likely to stay in the neighbourhood of theslowly varying equilibrium solution as long as a > cµ2/3 where c is someconstant.For α ≥ 43 and β ≥ 2, there’s no choice of γ that yields non-trivialapproximations, which means (3.5) may not be a good approximation. InSection 3.3, we find the approximation of the PDF by patching the WKBapproximation with the solution of a convection-diffusion process for α ≥ 43and β ≥ 2.3.2 WKB ApproximationTo focus on the early tipping case with larger amplitude of the noise relativeto the square root of the drifting rate (µ < ǫ2, same as β > 2), we need to513.2. WKB Approximationconsider two different range of bifurcation parameter a, as a = ǫαs with α >43 or α < 43 . Here we focus on finding the PDF for β ≥ 2 and α < 43 wherethe leading order approximation of the PDF is dominated by the potentialwell of the system (3.1). In Section 3.3, we discuss the approximation of thePDF of the system (3.1) given by the WKB approximation and a convection-diffusion process for α ≥ 43 and β ≥ 2. Table 3.1 provides the summary ofthe methods and scaled parameters we used based on different relationshipsof the drifting rate, noise, and bifurcation parameter. Table 3.2 summarizesthe parameters and variables we used in Section 3.2 and Section 3.3 as areference.µ a xWKB O(ǫβ) O(ǫα), α < 23β, α < 43 O(ǫα/2)Gaussian O(ǫβ), β < 2 O(ǫα), α > 23β O(ǫα/2)Convection-diffusion O(ǫβ), β > 2 O(ǫ4/3) O(ǫ)Table 3.1: Summary of the scaling parameters and corresponding methodsRecalling the scaling for (3.7), we substitute γ = 2 − 32α into (3.6) andget− ψs = ǫ2−β · 2ξ − ǫ32α−β(s− ξ2)ψξ + ǫ32α−β · 12ψξ2 + ǫ2−β · 12ψξξ. (3.11)By substituting the general expansion for ψ(ξ, s) in (3.11) asψ(ξ, s) ∼ ψ0 + ǫδ1ψ1 + ǫδ2ψ2 + · · ·, (3.12)where δj , j = 1, 2, 3 · · · are determined as shown below, we have the leading523.2. WKB Approximationorder contribution ψ0 satisfies (3.7) at O(ǫ32α−β)ψ0 = 2sξ −23ξ3 + f0(s), (3.13)where f0(s) is an undetermined function of time.We chooseδj = j · (2−32α), j = 1, 2, 3, · · · , (3.14)and β + 32α < 4 in order to get the following form in each order beyond theleading orderO(ǫ 32α−β+δj ) : (s− ξ2)ψjξ − ψ0ξψjξ = Gj(ψ0, ψ1, ..., ψj−1), j = 1, 2, 3, · · · ,(3.15)where ψ0, ψ1,..., ψj−1 are known from previous balancing procedures. Bysubstituting (3.13) into (3.15) we have(ξ2 − s)ψjξ = Gj(ψ0, ψ1, ..., ψj−1), j = 1, 2, 3, · · · . (3.16)One of the singularities in (3.16) is ξ = √s which conflicts with x = √awhich is the stable equilibrium of (3.1) for a a constant. To eliminate thesingularity so that the solution does not grow to ∞, we apply the secularitycondition asGj(ψ0, ψ1, ..., ψj−1)|ξ=√s = 0, j = 1, 2, 3, · · · . (3.17)533.2. WKB ApproximationFor j = 1, 2, we haveO(ǫ 32α−β+δ1) = O(ǫ2−β) :(ξ2 − s)ψ1ξ = G1 = ǫβ−2ψ0s + (2ξ +12ψ0ξξ), (3.18)O(ǫ 32α−β+δ2) = O(ǫ4−β− 32α) :(ξ2 − s)ψ2ξ = G2 = ǫβ−2ψ1s +12ψ1ξξ +12(ψ1ξ)2. (3.19)Terms ǫβ−2ψ(j−1)s, j = 1, 2 in (3.18) and (3.19) are correction terms inO(ǫ 32α−β+δj ) because its order is lower than O(ǫ 32α−β+δj+1), which makes usavoid extra expansions other than (3.14) and (3.15).By applying the secularity condition (3.17) on (3.18), we get(2ξ + f0s)|ξ=√s = 0, ⇒ f0 = −43s3/2 + C0, (3.20)ψ1ξ = ǫβ−2−(2ξ + f0s)s− ξ2 , ⇒ ψ1 = ǫβ−2 · 2ln(√s+ ξ) + f1(s), (3.21)where C0 is a constant and f1(s) is a function of s and needs to be deter-mined.By substituting (3.13), (3.20), and (3.21) into (3.19) and applying thesecularity condition (3.17), we havef1 =14lns− ǫβ−2lns+ C1, (3.22)where C1 is a constant that is related to the normalization of the PDF. The543.2. WKB Approximationapproximation of PDF p(ξ, s) is thenp(ξ, s) ∼ exp{ǫ 32α−2(2sξ − 23ξ3 − 43s3/2 + C0)+ ǫβ−2 ·[2ln(√s+ ξ)− lns]+ 14lns+C1 + · · · }.We change variables back by (3.3) and take C = exp{ǫ 32α−2 · C0 + C1}toget the approximation for p(x, t)p(x, t) ∼ C · exp{ 1ǫ2 (2ax−23x3 − 43a3/2) + 14ln(a)+ ǫβ−2[2ln(√a+ x)− ln(a)]}, (3.23)where a = a0 − µt.To estimate C, for a > 0 and far from its critical value 0, p(x, t) is wellapproximated by a normal distribution with mean √a. Introducing a Taylorexpansion with center x = √a in the exponential term of (3.23), we havep(x, a) ∼ C · exp(4a3/23ǫ2)· exp[−4√a2ǫ2 (x−√a)2 +O(ǫ)], (3.24)which implies C = 2√2πǫ2 · 4−ǫβ−2 for normalization.Figure 3.2 shows the PDF given by (3.23), which agrees with the resultsbased on numerical simulations given by the Euler-Maruyama method, com-pared with the PDF given by Gaussian approximation (3.10). For a = 0.05553.2. WKB Approximation−0.1 0 0.1 0.2 0.3 0.4 0.5 0.600.511.522.533.544.5xp(x,a)ǫ= 0.1; β= 3; µ= 0.001; a= 0.1Figure 3.2: Black dashed curve is the asymptotic approximation of the PDFgiven by (3.23). Magenta dash-dotted curve is based on the Gaussian typeapproximation (3.10). Red circle represents the numerical histogram under2.5 × 106 realizations by the Euler-Maruyama method.563.2. WKB ApproximationDescriptionx System variable in canonical model (3.1)a Slowly varying bifurcation parameter in canonical model (3.1)µ Drifting rate of a in canonical model (3.1)ǫ Amplitude of the white noise in canonical model (3.1)p(x, t) PDF of (3.1) in FPE(3.2)s Rescaled time variable as a = ǫαs for WKB approximationξ Rescaled system variable as x = ǫα/2ξ for WKB approximationα Defined in a = ǫαs to measure the distance between a and 0Defined in µ = ǫβ to capture the relationship betweenβ the drifting rate and the amplitude of the noise, in order toidentify different behaviours of (3.1)ψ(ξ, s) Defined in WKB form of the solution p(ξ, s) ∼ exp{ψ(ξ, s)/ǫγ}γ Defined in WKB form of the solution p(ξ, s) ∼ exp{ψ(ξ, s)/ǫγ}to identify different behaviours of (3.1)F (t) Solution of (3.1) with ǫ = 0 for Gaussian approximationg(t) Function related to F (t) for Gaussian approximationgn(t) Function related to F (t) for Gaussian approximationδj j = 1, 2, 3 · · · , exponents for ǫ in general expansion of ψfj(s) j = 0, 1, 2 · · · , function of time as part of the approximationof ψ(ξ, s)Gj j = 1, 2, 3 · · · , RHS of each order for the secularity conditionC Normalization coefficient for the PDFaˆ Value of a of O(ǫ4/3) for the tail approximationη Rescaled system variable as x = ǫη for the tail approximationxs Used to account for the diffusion process for x < xsxm Used to combine the WKB approximation and the tailapproximationxd The tail approximation is no longer valid for x < xdTable 3.2: Summary of the Parameters and Variables573.2. WKB Approximationrelatively far from its critical value, the heavy left tail of the PDF indicatesthat the early tipping is likely to happen while a keeps decreasing beforearriving the critical value.Summary of results: The method of multiple scales and WKB approx-imations shows the different leading order contributions of the FPE corre-sponding to the relationship of the drifting rate of the bifurcation parameterand the amplitude of the noise, which reflect different behaviours of the sys-tem (3.1). When β < 2, which is the same as µ > ǫ2, the Gaussian-typeapproximation given by [34] is applicable, which indicates the sample pathsof the system is likely to follow the slowly varying equilibrium solution ofthe system and no early tipping happens. When β > 2, which is the same asµ < ǫ2, the leading order terms of the WKB approximation is driven by thepotential well of the system at every value of the bifurcation parameter. Thedifferent behaviours of the system (3.1) separated by the relationship µ = ǫ2is consistent with the conclusion of [28] based on a sample-path analysis: ifnosie ǫ is larger than the square of the drifting rate of the bifurcation pa-rameter µ1/2, the sample paths are likely to escape from the neighbourhoodof the stable branch for positive a of order ǫ4/3; if the amplitude of the noiseǫ is smaller than the square root of the drifting rate µ1/2, the sample pathsare likely to stay in the neighbourhood of the slowly varying equilibriumsolution as long as a > −c · v2/3 where c is some constant. However, ourapproach provides an explicit expression for the approximation of the PDF,which has not been given in previous studies [28, 38].583.3. PDF Tail Approximation3.3 PDF Tail Approximation for TrajectoriesCrossing the Unstable BranchThere are two constraints for the WKB approximation (3.23): (1) β+ 32α <4; (2) x > −√a. For β + 32α ≥ 4 which indicates β − 2 ≥ 2 − 32α, theterm of O(ǫβ−2) in (3.18) and (3.19) should be moved to the higher ordercontributions as ǫβ−2 ≤ ǫ2− 32α. For x = −√a, (3.23) becomesp(x, t) ∼ C ·(√a+x)2ǫ−4exp{ 1ǫ2 (2ax−23x3−43a3/2)+14ln(a)−ǫβ−2ln(a)} = 0.For a closer to the bifurcation point where the stable and unstable equi-librium are close to each other, the approximation (3.23) is not valid forx < −√a since it exponentially increases while x decreases from −√a, whichconflicts with the results from the fact that the PDF approaches 0 while xapproaches −∞.For a closer to the bifurcation point, there is a non-trivial probabilityfor the system to remain near the unstable branch, as there is a slowlyvarying equilibrium solution for the deterministic system (3.1) (ǫ = 0) nearthe unstable branch, described in [28, 32]. The probability for the systemstaying in the neighbourhood of the unstable equilibrium cannot be ignoredfor estimating the probability of early tipping which is the probability forx far below the unstable equilibrium −√a. Hence, it’s important to get anidea about what happens near x = −√a.For a = aˆ = O(ǫ4/3) where the approximation (3.23) is still valid since593.3. PDF Tail Approximationα > 43 , we divide the range of x into two parts by introducing another pa-rameter xm, which is close to but larger than −√aˆ. For x > xm, (3.23)provides a reasonable approximation as x is attracted by the stable equi-librium√aˆ when it is in the neighbourhood of the stable equilibrium. Forx ≤ xm, we determine the local behaviour of (3.2) using a scaled variablex = ǫη, yielding the equation for p(η, t)pt = −aˆǫ pη + ǫ(η2pη + 2ηp) + ǫβ−1tpη +12pηη , (3.25)p(η, 0) = pqs(η, aˆ), η < xs/ǫwhere aˆ = O(ǫ4/3) is where the WKB approximation (3.23) is still valid, andpqs is the asymptotic approximation (3.23). We choose xs > xm in order toaccount for the diffusion process given by (3.25) that influences x ≤ xm.By substituting the general expansion for p(η, t) in (3.25) asp(η, t) ∼ p0 + ǫp1 + ǫ2p2 + · · · , (3.26)we find the leading order equation for p0O(1) : p0t +aˆǫ p0η =12p0ηη , (3.27)p0(η, 0) = pqs(η, aˆ),which is a convection-diffusion equation with no boundary conditions [33].603.3. PDF Tail ApproximationIts solution based on original variables x and t is given byp0(x, t) =∫ ∞−∞1√2πtexp[−(x− aˆt− z)22t]pqs(z, aˆ)dz. (3.28)To obtain higher order contributions, we proceed similarly to the analysis(3.14) - (3.15). We assume ǫ2 < aˆ < ǫ and 2 < β < 3, which covers thecases β > 2, α < 43 and β > 2, α > 43 , where we cannot use the WKBapproximation. We get the following general form for the equations of thehigher order corrections.O(ǫj) : pjt +aˆǫ pjη −12pjηη = Fj(pj−1, η, t), j = 1, 2, 3, · · · (3.29)Fj(pj−1, η, t) = η2p(j−1)η + 2ηpj−1 + ǫβ−2tp(j−1)η,where pj−1 is known from lower order equations of O(ǫj−1), j = 1, 2, · · · .Writing in terms of S = t and x = η − aˆǫ t, we have a non-homogenous heatequation in O(ǫj) with zero initial conditions and the solution is given bypj(x, t) =∫ t0∫ ∞−∞1ǫ√2π(t− s)exp[−(x− aˆt− z)22ǫ2(t− s)]Fj(pj−1(z, s), z, s)dzds(3.30)with Fj(pj−1(x, t), x, t) = 1ǫ (x2p(j−1)x + 2xpj−1) + ǫβ−2tp(j−1)x.The approximation (3.28) and (3.30) is valid for x close to 0. We intro-duce a new parameter xd < −√aˆ which is relatively far from 0. For x < xd,613.3. PDF Tail Approximation(3.28) and (3.30) is no longer valid since for these values of x the trajectoryis in the range where the convection dynamics dominate over the diffusion,which drives the system strongly away from the unstabie equilibrium. Thisbehaviour is equivalent to the trajectory of a particle that is well below theunstable equilibrium, with negligible probability of returning to it. In fact,if we focus on x < 0 and |x| = O(1) in the FPE (3.2), the leading ordercontributions are p0t − x2p0x = 2xp0, which is a quasi-linear PDE, of whichthe solution propagates to −∞ along time. Here we choose xd = −√aˆ− 2ǫwhere the approximation given by (3.28) and (3.30) exponentially increasesinstead of decreases to 0 for decreasing x.To get the PDF for x > xd at a = aˆ − µt for a certain value of t, wefind the approximation for xd < x < xm and x > xm separately and patchthem by finding a constant determined by normalization. For xd < x < xm,we use (3.30) on xd < x < xs with initial condition given by (3.23) ata = aˆ to get the partial PDF and the probability for xd < x < xm asP (xd < x < xm) =∫ xmxd p(x, a)dx. For x > xm, the PDF is approximatedby (3.23). The normalization coefficient C in (3.23) satisfies that P (x >xm) =∫∞xm p(x, a)dx and P (x > xm) + P (x < xm) = 1. Here P (x < xm) =P (xd < x < xm)+P (x < xd) and the probability P (x < xd) =∫ xd−∞ p(x, t)dxis given by the PDF (3.28) and (3.30), since below xd the system (3.1) isconsidered “escape” based on the above analysis. This approximation is notvalid for large values of t as the WKB approximation (3.23) moves alongwith the stable equilibrium while a slowly decreases and the PDF controlledby the convection-diffusion process increases. The reason that we can patchthe solutions (3.23), (3.28) and (3.30) together is that they are the solutions623.3. PDF Tail Approximationof the same FPE (3.2). Even though they are of different orders of variables,they overlap for a certain range of x.To approximate the PDF for larger t, for example for a = aˆ − 2µt, weinitiate the same approximation and patching procedure based on the PDFat a = aˆ−µt. The related xd, xs and xm depend on ǫ, β and t. The iterationcan be described as follows:(1) Given aˆj = aˆ−j ·µt, and corresponding xdj , xmj and xsj, j = 1, 2 · · · ,we can find the PDF of (3.1) p(x, aˆj) based on p(x, aˆj−1). How xdj , xmj andxsj are chosen is described below. For j = 1, the approximation for xd1 <x < xm1 is given by (3.28) and (3.30) with initial condition at a = aˆ andthe approximation for x > xm1 is given by (3.23) in which the normalizationcoefficient C is determined by the probability P (x < xm1)+P (x > xm1) = 1.(2) For j = 2, 3, · · · , the approximation for xdj < x < xmj is givenby(3.28) and (3.30) with initial condition at a = aˆj−1 for x < xsj. Theapproximation for x > xmj is given by (3.23).(3) To determine the normalization coefficient in (3.23), we use P (x <xmj) + P (x > xmj) = 1, where P (x > xmj) =∫∞xmj p(x, aˆj)dx. UnlikeP (x < xm1) in (1), P (x < xmj) =∫ xmjxdj p(x, aˆj)dx +∑j−1k=1 P (x < xdk), inorder to capture the probability of early tipping.Hence, even though we don’t have an explicit expression of the PDF forall the values of x, we are able to estimate the probability of early tipping ataˆj , j = 1, 2, 3 . . . by estimating the probability P (x > −√aˆj), where −√aˆjrepresents the unstable equilibrium at a = aˆj .In order to get a better approximation, we need to choose xsj > xmj >−√aˆj > xdj and t properly. To consider the effect of the convection-diffusion633.3. PDF Tail Approximationprocess thoroughly, we need a reasonable range of x in the initial conditionto describe the evolution of the process locally, which gives the evolution ofthe probability density, for trajectories near xmj , including x < xmj . Thelarger xsj is, the larger xmj can be for the patching between the left tailof the PDF and the WKB approximation. To balance the accuracy andthe expense of the computation time, we choose xsj = −√aˆj + 10ǫ in theapproximation. Similar to the choice of xd, we choose xdj = −√aˆj − 2ǫ,below which we assume x escapes. Since t represents the time length of eachiteration step, it also affects the value of xmj . Based on the exponents in(3.28), if t is too small, the effect of the convection-diffusion process is alsosmall. If t is too large, the initial PDF under the effect of the convection-diffusion process moves to the far left which can not be patched with theWKB approximation. By considering the value of exp(− aˆ2j−1t22t ) in (3.28) ineach iteration, we choose t = 2 and xmj = −√aˆj + ǫ for the approximation.Figure 3.3 shows how the iteration works to provide the approximation ofthe PDF for (3.1). The values of xd, xm and xs move along with the locationof the unstable equilibrium. For a closer to the critical value ae = 0, theprobability of x < −√a isn’t negligible as the approximation of the earlytipping.Figure 3.4 shows the PDF given by (3.23) and (3.30) agrees with theresults from numerical simulations given by the Euler-Maruyama method forx > xd. But the PDF given by (3.10) severely underestimates the probabilityfor x below the unstable equilibrium.Summary of results: For smaller a where the stable and unstable equi-librium are both closer to 0 and no longer far from each other, the WKB643.3. PDF Tail Approximation−0.4 −0.2 0 0.2 0.401234567a= 0.031056xp(x,a)−0.4 −0.2 0 0.2 0.401234567a= 0.020993xp(x,a)−0.4 −0.2 0 0.2 0.401234567a= 0.010931xp(x,a)−0.4 −0.2 0 0.2 0.401234567a= 0.0053409xp(x,a)Figure 3.3: The PDF approximation for (3.1) with ǫ = 0.05, β = 2.5,given by (3.23), (3.28) and (3.30) at different a(t). The vertial blue dottedline represents the unstable equilibrium −√a. The solid red circle/bluetriangle/black circle represents the corresponding xd/xm/xs. xs > 0.4 inLower right figure, which is shown in Figure 3.4.653.3. PDF Tail Approximation−0.4 −0.2 0 0.2 0.4012345xp(x,a)a= 0.0053409Figure 3.4: PDF of (3.1) for ǫ = 0.05, β = 2.5 at a = 0.005 where tippingis likely to happen. Black dashed curve is the asymptotic approximation ofthe PDF given by (3.23), (3.28) and (3.30). Magenta dash-dotted curve isbased on the Gaussian type approximation (3.10). Red circle represents thenumerical histogram under 2.5 × 106 realizations by the Euler-Maruyamamethod. Blue dotted vertical line represents the unstable equilibrium for(3.1). The solid red circle/blue triangle/black circle represents the corre-sponding xd/xm/xs.663.3. PDF Tail Approximationapproximation is not valid near the unstable equilibrium. A local analysisof the system (3.2) near x = 0 and a = 0 is applied and the system turnsout to be a series of convection-diffusion equations. By finding the solu-tion asymptotically and combining it with the WKB approximation nearthe stable equilibrium through normalization, we find the approximationof the PDF of the system (3.1) which provides the probability of the earlytipping for any a > 0. For a approaching 0, nonzero PDF rises near theunstable equilibrium and the probability of the system below the unstableequilibrium is no longer negligible, which indicates that it is possible for thesystem to escape from the stable state for those a > 0. By the combinationof the WKB approximation and the convection-diffusion process near theunstable equilibrium, we provide the approximation of the PDF under theinfluences of both the drifting bifurcation parameter and the additive noisefor tipping before the system reaching the bifurcation point, which is notgiven in previous studies [17, 28, 38].67Chapter 4Applications4.1 Morris-Lecar ModelWe consider the Morris-Lecar(ML) model [10, 45] with a slowly increasinginput current and an additive periodic external input. The ML model hasbeen used widely as a two-dimensional model, relatively straightforward toanalyze, which captures certain dynamical behaviour observed for differenttypes of neurons with a variety of states. It has the general formγdvdtˆ= −gCamˆ∞(v)(v − vCa)− gK(v − vK)w − gL(v − vL) + Ibias + Iext,(4.1)dwdtˆ= κˆ(v)(wˆ∞(v)− w),dIbiasdtˆ= µˆ,mˆ∞(v) =12(1 + tanh(v − v1v2)), wˆ∞(v) =12(1 + tanh((v − v3v4)),κˆ(v) = φcosh(v − v32v4),where v represents the membrane voltage potential and w is the probabilityof opening of K+ channels (0 ≤ w ≤ 1). The periodic external input isIext, and Ibias represents the bias or base current in the neuron. We take684.1. Morris-Lecar Modelthe commonly used fixed constants as given in [45]: vK = −84mV, vL =−60mV, vCa = 120mV, c = 20µF/cm2, v1 = −1.2mV and v2 = 18mV,gK = 8mS/cm2, gL = 2mS/cm2, gCa = 4.4mS/cm2, v3 = 12mV, v4 =17.4mV, and φ = 115ms−1. The dynamical behavior of (4.1) has been wellstudied for Iext = 0, with the current Ibias = constant often playing therole of the bifurcation parameter. For that case there is a SNIC (saddle-node on an invariant circle) bifurcation [44] at Ibias = Ic ≈ 44.09µA/cm2,and corresponding voltage v = vc ≈ −27.14mV, so that for Ibias < Ic.there exist a stable equilibrium and unstable equilibria. For Ibias > Ic,the system exhibits oscillations whose frequency increases from zero withIbias − Ic, and there is a Hopf bifurcation for Ibias above 60µA/cm2, whichwe do not consider here. The bifurcation structure of a Morris-Lecar neuronunder a sinusoidal stimulus has been studied [36, 37], where the bifurcationparameter Ibias is a constant. Different types of bifurcations occur at thephase boundary of the non-firing state, such as the inverse flip bifurcation,the periodi-doubling bifurcation and the saddle-node bifurcation.We focus on transitions from the stable steady state to oscillations fora current that is slowly increasing with rate µˆ≪ 1, modeled here as slowlyvarying Ibias for simplicity but alternatively it could be an input current, anda periodic forcing Iext = Aˆ sin(Ωˆt). Sample trajectories superimposed on thebifurcation structure are shown in Figure 4.1. To determine the location ofthe tipping points, that is, value of Ibias for the transitions, we transform(4.1) to a form similar to the canonical model (2.1), rescaling time tˆ = γtand normalizing the system about the bifurcation point (vc, Ic),694.1. Morris-Lecar Modeldxdt = h(x)− gK(x+D)w + b(t) +Asin(Ωt)= b(t)− gCam∞(x)(x+ 1−vCavc)− gL(x+ 1−vLvc) + Icvc,− gK(x+ 1−vKvc)w +Asin(Ωt), (4.2)dwdt = κ(x)(w∞(x)− w),dbdt = −µ = −γ|vc|µˆ,x = v − vcvc, b = Ibias − Icvc, Ω = γΩˆc,vcA = Aˆ, m∞(x) = mˆ∞(vc(x+ 1)).κ(x) = γκˆ(vc(x+ 1)), w∞(x) = wˆ∞(vc(x+ 1)), (4.3)For µ = 0 and A = 0, the bifurcation point is at the origin for the trans-formed system (4.2). Figure 4.1 shows the bifurcation diagram for (4.1) and(4.2), with trajectories for A = 0 and A 6= 0 obtained from the correspond-ing systems with µ 6= 0. We note here that for the transformed system,the bifurcation parameter b and dependent variable x take O(1) values. Forthe system in this form, it is then straightforward to identify different scalesrelated to slow drift, high and low frequency, or large amplitude. Such anidentification is important in choosing the correct approach to approximatethe tipping point.Before considering the case with A 6= 0, we give an analytical expressionfor the tipping point without oscillatory forcing A = 0, which is comparedbelow with the tipping point for A 6= 0. It is useful to find the tipping point704.1. Morris-Lecar Model0 10 20 30 40 50 60−60−40−2002040Ibiasv−0.5 0 0.5 1 1.5−2.5−2−1.5−1−0.500.51bxFigure 4.1: LEFT: The bifurcation diagrams of (4.1) indicating the stable(red solid) and unstable (red dashed) branches. The trajectories for Aˆ = 0(green dash-dotted line) and Aˆ 6= 0 (blue dashed line) are superimposedon these curves , with the corresponding tipping points indicated by solidcircles. Note that tipping occurs for a smaller value of Ibias for the caseAˆ 6= 0. RIGHT: Same as LEFT, but for the transformed system (4.2).714.1. Morris-Lecar Modelfor the generic system of the formdxdt = Da+ k0 + k1x+ k2x2 , dadt = −µ. (4.4)Applying the scaling of [32] described in Section 2.1, completing the squarein (4.4) and using a straightforward transformation again yields an Airyequation, from which we obtain the singularity corresponding to the tippingpoint for (4.4)ag = (D|k2|)−1/3ad −asD for as = k0 +k214|k2|, (4.5)where ad is the tipping point in (2.4) for the symmetric case and as givesthe relative shift for a general quadratic polynomial in (4.4) if D = 1.For A = 0 in (4.2) we obtain an approximation to the tipping point usinga local expansion about x = 0. Noting that 0 < κ(0) = O(1) we expect wto relax to w∞, so that the local approximation for (4.2) isdxdt ≈ b+ h00 + h01x+ h02x2− gK(x+D)[w∞(0) + w′∞(0)x + w′′∞(0)x2/2], (4.6)dbdt = −µ, h0j = j!h(j)(x)|x=0.724.1. Morris-Lecar ModelThen the approximate tipping point for A = 0 is given by (4.5) withk0 = h00 − gKDw∞(0),k1 = h01 − gK(w∞(0) +Dw′∞(0)), (4.7)k2 = h02 − gK(w′∞(0) +Dw′′∞(0)2).Note that k0 = k1 = 0 in (4.7) following from the quadratic form for thesaddle-node bifurcation at b = 0 and x = 0, so that resulting tipping pointfor (4.6) is bml = |k2|−1/3ad.4.1.1 High frequency oscillation Ω ≫ 1We use the approach of Section 2.2.2 to study the case µ≪ 1 and Ω ≫ 1 in(4.2), again writing Ω = µ−λ. We restrict our attention to the case whereA = O(1).We introduce the fast time scale T = µ−λt and the slow time scaleτ = µt, and substitute the multiple scale approximations for x and wx ∼ x0(T, τ) + µλx1(T, τ) + · · · , (4.8)w ∼ w0(T, τ) + µλw1(T, τ) + · · · , (4.9)into (4.2), yielding the sequence of equationsO(1) : x0T = w0T = 0,⇒ x0 = x0(τ), w0 = w0(τ), (4.10)O(µλ) : x1T = b+ h(x0)− gk(x0 +D)w0 +Asin(T ),w1T = κ(x0)(w∞(x0)− w0). (4.11)734.1. Morris-Lecar ModelWe apply a solvability condition similar to (2.9), namely that the innerproduct of the right-hand side of (4.11) with the homogeneous solutions (xand w constant with respect to T ) vanish. This yieldslimL→∞1L∫ L0[b+ h(x0)− gK(x0 +D)w0 +A sin(T )] dT = 0, (4.12)limL→∞1L∫ L0κ(x0)(w∞(x0)− w0) dT = 0.Then the asymptotic approximation of x for b(τ) = O(1) isx ∼ x0 + µλ[−Acos(T )] + · · · , where − b = h(x0)− gK(x0 +D)w∞(x0)(4.13)w ∼ w∞(x0) + · · · ,where x0 in (4.13) is the stable equilibrium of (4.2), shown as the solidcurve in Figure 4.1 (RIGHT). The asymptotic approximation describes theattracting solution of (4.2) away from b = 0, but does not give the tippingpoint.To determine the location of the tipping point, we use a local approxi-mation for x≪ 1 similar to that used in Section 2.2.2,x = X (T ) + z(T, s), X (T ) = −Ω−1Acos(T ), (4.14)w = W0(T, s) + Ω−1W1(T, s),z = z0(T, s) + Ω−1z1(T, s), s = µ1/3t, (4.15)that is, the solution is given by oscillations on the T scale plus a correction744.1. Morris-Lecar Modelz. We can introduce additional scalings of z = O(µ1/3) and b = O(µ2/3) asin (2.2.2), but as it does not affect the result, we use z and b for simplicity.We also introduce quadratic polynomials given by Taylor expansions of h(x),κ(x) and w∞(x) about x = 0, facilitating explicit expressions for the tippingpoint. Then substituting (4.14)-(4.15) in (4.2), we find that W0T and z0Tvanish, so that W0 =W0(s), z0 = z0(s). Approximations for W0 and z0 areobtained at the next order, analogous to (2.20) as shown in Appendix D.We apply the solvability condition as in (4.12), to the equations for W1 andz1 yielding for W0W0 ∼ w∞(0) +W00 +W01z0 +W02z20 ,with the coefficients W00, W01, and W02 given in (D.3). Then the equationfor z0 written in terms of t is,z0t = b+ h00 + h01z0 + h02z20 + h02A22Ω2 − gK(z0 +D)W0. (4.16)Applying (4.4)-(4.5) yields the singularity for (4.16)bhf = |hhf2 |−1/3ad − bs for bs = hhf0 +(hhf1 )24|hhf2 |, (4.17)hhf0 = −gkDW00 + h02A22Ω2 , hhf1 = −gK(W00 +D(W01 − w′∞(0))),hhf2 = h02 − gK(W01 +DW02),where we have used the fact that k0 and k1 in (4.7) vanish. Note that forΩ ≫ 1, W01 − w′∞(0) and W00 are proportional to Ω−2, as are hhf0 and hhf1 .754.1. Morris-Lecar ModelThen we recover (4.7) as Ω−2 → 0.Ibias20 25 30 35 40 45 50λ0.30.40.50.60.70.80.911.11.2Ibias20 25 30 35 40 45 50λ0.30.40.50.60.70.80.911.11.2Figure 4.2: The value Ibias corresponding to tipping in (4.1) vs. λ for Ω =µ−λ. The asymptotic results obtained from (4.17) (solid lines) are comparedwith numerical simulations (dash-dotted lines), obtained using the conditionfor tipping as v = 27.13 in (4.1) by the 4th order Runge-Kutta method.LEFT: Results are shown for amplitude A = 2 (Aˆ = −54.275) and twodrifting rates µ = 0.001 (µˆ = 0.0007) (red) and µ = 0.02 (µˆ = 0.0147)(blue) with asymptotes for large Ω approximated by (4.7) (vertical dottedlines) as Ibias = 47.65 and Ibias = 44.58, respectively. RIGHT: Results areshown for µ = 0.003 (µˆ = .0022) and amplitude values A = 2 (red) andA = 8 (Aˆ = 217.10) (blue), with earlier tipping for larger ratios A/Ω.In Figure 4.2 we compare the asymptotic approximation for the tippingpoint for Ibias obtained from (4.17) as a function of the exponent λ forΩ = µ−λ with numerical simulations by the 4th order Runge-Kutta methodfor different drifting rates µ and amplitudes A. For larger values of theratio A/Ω we see substantial advance of the tipping point. For larger λcorresponding to higher frequencies Ω, the delay effect caused by slowlyvarying drifting rate µˆ dominates the location of the tipping point, regardlessof the amplitude of the oscillatory forcing. For smaller λ, larger amplitudeAˆ triggers ealier tipping.764.1. Morris-Lecar Model4.1.2 Low frequency oscillation Ω ≪ 1As in Section 2.3, we write Ω = Cµ, rescale time as τ = µt and write x andw as functions of b(τ) = b0 − τ . Then (4.2) becomes−µdxdb = b+ h(x)− gK(x+D)w +Asin(C(b0 − b)), (4.18)−µdwdb = κ(x)(w∞(x)− w). (4.19)Substituting the asymptotic expansion (4.9), we find the leading order ap-proximation of w(x) is w0 = w∞(x0) for x0 ∼ O(1) and for x0b+ h(x0)− gK(x0 +D)w∞(x0) +Asin(C(b0 − b)) ≡ Φ(x0, b) = 0. (4.20)As in Section 2.3, we expect the tipping to occur in the region where theouter approximation is not valid, for a pair (br, xr). A Taylor expansion ofΦ about xr evaluated at br, keeping up to quadratic terms, yieldsx0 − xr = −12Φr2[Φr1 ±√(Φr1)2 − 4Φr2 · F(br)], Φrj =∂jΦ∂xj (xr, br),F(b) = Φr0 − br −Asin(C(b0 − br)) + b+Asin(C(b0 − b)), (4.21)where the positive sign corresponds to the solution that is attracting. Itfollows that the approximation for x0 is no longer valid for b < br where4F(br) = (Φr1)2/Φr2 and 2(x0 − xr) = −Φr1/Φr2. As x0 → xr, Φr1 → 0, fromwhich we conclude that xr = xc = 0, the value of x at the saddle-nodebifurcation point for the static case. As a result we expect that tippingoccurs near b = br and x = 0. Note that F(b) then has the same functional774.1. Morris-Lecar Modelform as f(a) in (2.26).As in Section 2.3 we use a local analysis in terms of b− br = µ2/3B andx = µ1/3χ into (4.18), yielding− µ2/3χB = F ′(br)µ2/3B + k2(µ1/3χ)2 + o(µ2/3), (4.22)where F ′(br) = 1− CA cos(C(b0 − br)), and k2 is the same as in (4.7). Thenusing the result (4.5) yields the tipping pointblf = br +ad(F ′(br)|k2|)1/3. (4.23)The result for blf relies on F ′(br) = O(1), and as in Section 2.3, if F ′(br)vanishes the approximation (4.23) is not valid. Then we again expect to finda jump in the location of the tipping point as a function of the amplitudefor a critical value A = Ac. Following Section 2.3 we provide a local analysisfor (b,A) close to (bm, Am), which satisfy the conditionsF(bm) = bm +Asin(C(b0 − bm)) = 0, (4.24)F ′(bm) = 1− CAmcos(C(b0 − bm)) = 0, ⇒ 1 + C2b2m = C2A2m. (4.25)To find Ac > Am for which there is a jump in the location of the tippingpoint, we consider the concavity of x in (4.18) with respect to b near bm.Substituting local variablesb− bm = µη, x = µξ, and A = Am + µA1, (4.26)784.1. Morris-Lecar Modelinto (4.18) and (4.19), we find the leading contributions to x′(b) and x′′(b)in terms of the local variables,d2ξdη2 = µA1CAm− µ(C2ηbm)− µk2ξ2A1bmAm+O(µ2), (4.27)dξdη =bmA1Am+O(µ), (4.28)for k2 from (4.7). From (4.27) we see that the line corresponding to theconcavity threshold, across which the trajectory of x changes its concavity,is given byξMLc = −C2Am2k2A1η + 12Ck2bm. (4.29)We compare the slope of (4.29) to the slope of x from (4.28), as in (2.35)-(2.37). Where these two slopes are equal, we find the change in concavity ofthe trajectory that corresponds to tipping. Setting the slope of (4.29) equalto (4.28), and solving for A1 yieldsA1 =CAm√2|k2|bm⇒ Ac = Am +ΩAm√2|k2|bm. (4.30)Figure 4.3(LEFT) compares the asymptotic approximation of the lo-cation of the tipping point given by (4.23) to the results from numericalsimulations by the 4th order Runge-Kutta method. We show the locationof the tipping point as a function of the amplitude Aˆ, with the jump inthis curve occuring at the amplitude threshold Ac given by (4.30). Figure4.3(RIGHT) illustrates the character of the early transition near A = Acfor (4.1) with early tipping occuring for A > Ac at but not for A < Ac.794.1. Morris-Lecar Model20 30 40 50 60050100150200ItipAˆ0 10 20 30 40 50 60 70−60−40−2002040IbiasvFigure 4.3: µ = 0.04, Ω = 0.2, I0 = −55. LEFT: Numerical simulations ofthe location of the tipping points of (4.1) (Red dashed line with dots) fordifferent Amplitude Aˆ based on the 4th order Runge-Kutta method, com-pared with the asymptotic approximation given by (4.23) (Red solid), withtipping defined as v = 27.14. The black circle shows the amplitude thresh-old Ac = 14.11 and the Ibias value corresponding to bm = 0.41. RIGHT:Numerical simulations for different amplitude Aˆ: −Aˆ = 13.02 (black solid)and −Aˆ = 15.20(blue dashed).804.1. Morris-Lecar ModelFor the value of Ibias where early tipping occurs, there is not an attractingoscillatory solution nearby, so following the tipping event, the system doesnot remain in the oscillatory state. This behavior is in contrast to the sys-tem (2.1), for which the tipping is irreversible. For some combinations of b0and A, the tipping point for (4.2) associated with Ac is closer to the SNICbifurcation, in which case the tipping for A > Ac the system remains in theoscillatory state.Summary of results: Our method developed in Chapter 2 is applicable tothe higher dimensional system with a saddle-node bifurcation, slowly varyingbifurcation parameter and additive periodic forcing, for example, the Morris-Lecar model with a slowly increasing input current and an additive periodicexternal input current. It is important to normalize the system (4.1) withproper scaling and apply the approaches for high or low frequency case tothe normalized system. For the high frequency oscillation, the location ofthe tipping point is determined by the local bifurcation structure near thebifurcation point, the slowly drifting bifurcation parameter, the amplitudeand the frequency of the periodic forcing. It is necessary to carefully analysethe slowly varying variable w in both outer and inner approximations to havea reasonable approximation for x. For the low frequency case, tipping occursnear values of b where the trajectory is near x ∼ 0, which is determined bynot only the varying bifurcation parameter, but also the amplitude, thefrequency and the initial condition of the external oscillatory forcing. Thediscontinuity of the location of the tipping point related to the amplitude ofthe periodic forcing is captured by the concavity analysis similar to the onein Chapter 2.814.2. Energy Balance Sea Ice Model4.2 Energy Balance Sea Ice ModelWe apply the analysis to an energy balance model developed by Eisenmanand Wettlaufer [16], as mentioned in the Introduction. The model describessea ice dynamics in terms of the change of the energy per unit area, storedeither as latent heat in sea ice or as sensible heat in the ocean mixed layer,E ≡−Lihi E < 0cmlHmlTml E ≥ 0,(4.31)where Li is the latent heat of fusion for sea ice, hi is the thickness of sea ice,cml is the heat capacity of the ocean mixed layer, Hml is the depth of themixed layer, and Tml is the temperature difference from the freezing point0◦C. The reduced system, derived from the full heat conduction equationsfor sea ice thermodynamics of [1], expresses the energy dynamics in terms ofsensible and latent heat fluxes, outgoing longwave and incoming shortwaveradiative surface fluxes. We focus on the partially linearized version of themodel in [16] that neglects ice export which has the formdEdt = [1− α(E)]FS(t)− F0(t) + ∆F0 − FT (t)T (t, E) + FB . (4.32)Here FS(t) describes the shortwave radiation flux which is seasonally vary-ing, with [1 − α(E)]FS(t) describing surface shortwave radiation related toalbedo feedback,α(E) = αml + αi2 +αml − αi2 tanh(ELihα), (4.33)824.2. Energy Balance Sea Ice Modelwhere αml and αi are constants for albedo feedback for the mixed layer andsurface ice, respectively. The factor hα is an ice thickness range for smoothtransition from αi to αml. The temperature dependence of the surface flux isobtained via a linearization of the Stefan-Boltzmann equation, resulting inthe outgoing longwave radiation contributions F0(t)+FT (t)T (t, E), with theseasonal variation in F0(t) and FT (t) determined from combined modelingand observational data (see SI Appendix of [16] for details). The partiallinearization in [16] is given by taking the surface temperature T (t, E) asT (t, E) = EcmlHml. (4.34)The term FB describes the heat flux into the bottom of the sea ice or theocean mixed layer. The additional contribution ∆F0 captures a change inthe surface heat flux, representing warming in the model for ∆F0 > 0.In [16] the nonlinear threshold behavior of the energy was studied throughthe bifurcation structure of (4.32) with the seasonally varying terms F0(t),FS(t) and FT (t). Attracting periodic states for fixed ∆F0 were determinedby simulating the equation, allowing the system to reach the attracting state,recording the minimum and maximum value of E for the particular value of∆F0. For the parameters in [16], there are two branches of attracting statesfor E > 0 and E < 0 in terms of ∆F0, the bifurcation parameter describingthe change to the surface heat flux. These states are bistable for a rangeof ∆F0, with each branch losing stability to the other via a saddle-nodebifurcation. Figure 4.4 shows the branch corresponding to the stable steadystate of the averaged version of (4.32), that is, with the seasonally varying834.2. Energy Balance Sea Ice Modelterms F0(t), FS(t) and FT (t) replaced by their averages, and for E ≥ 0and ∆F0 a positive constant. We denote the saddle-node bifurcation valueon this branch as (∆F0c, Ec), and also show the unstable branch for thatsystem. As would be expected from the analysis of previous sections, thesaddle-node bifurcation for the averaged system differs considerably fromthe tipping point when there is a large amplitude periodic forcing.In considering the case where the surface heat flux varies slowly overtime with a rate µ˜ we use the asymptotic approaches described in previ-ous sections to study tipping from the attracting branch E ≥ 0, notingthat by symmetry a similar approach applies for the other branch E < 0.We contrast the location of the tipping point for the averaged model withthe tipping location in the case where the oscillatory functions F0(t), FS(t)and FT (t) are not replaced with their averages. Since we obtain analyticalexpressions for the tipping point, we can easily explore the parametric de-pendence on the location of the tipping point, which has implications forthe bifurcation structure as a function of the parameters. We note that thecomputational approach used in [16] to explore bifurcations corresponds tothe limit µ˜→ 0.The fluctuations in F0(t), FS(t) and FT (t) are given as a series of monthlymeasurements [16] (SI Appendix), which we interpolate and approximate us-ing a (finite) Fourier series approximation. We assume that the fluctuationsare the same year to year. To allow easy identification of the oscillatoryforcing terms, we write (4.32) in a form analogous to that of (2.1), shiftthe bifurcation point for the averaged system from (∆F0c, Ec) to the originas shown in Figure 4.4. We rescale the system appropriately, so that the844.2. Energy Balance Sea Ice Modelresulting dependent variable x and corresponding slowly varying parameterb are O(1) quantities. Substitutingx = E − EcEcb = ∆F0 −∆F0cEc(4.35)in (4.32), and taking ∆F0 to vary slowly in time with rate µ˜ yieldsdxdt = g1(t) + g2(t)tanh(g3(x+ 1)) + g4(t)x+ b, (4.36)dbdt = −µ = −µ˜/Ec.g1(t) = E−1c[(1− αml + αi2 )FS(t)− F0(t) + FB]+ ∆F0cEc+ g4(t),g2(t) = −αml − αi2EcFS(t), g3 =1LihαEc, g4(t) = −FT (t)cmlHml.In contrast to (2.1), which has only an additive periodic forcing, (4.36)has both multiplicative oscillatory forcing terms with coefficients g2(t) andg4(t) and an additive oscillatory term g1(t). So we consider approximationsfor the terms involving g2 and g4 that replace the multiplicative oscillationswith reasonable additive oscillations. The oscillations in g1 have an ampli-tude of ≈ 20, g2(t) has an amplitude of ≈ 5, while g4(t) has oscillations withrelatively small amplitude of ≈ 0.5. We anticipate that the fluctuations ofg4(t), as compared with those of g1(t) and g2(t), have a negligible effect ontipping, generally occurring for x = O(1) or smaller. So we neglect the oscil-lations in the term g4(t)x by replacing g4(t) with its average. Furthermorethe results from the previous sections suggest that large oscillations in g1drive tipping for O(1) values of the bifurcation parameter b and thus forO(1) values of x, for which tanh(g3(x + 1)) ≈ 1. Then, it is reasonable to854.2. Energy Balance Sea Ice Modelapproximate g2(t)tanh(g3(x+1)) ≈ G2tanh(g3(x+1)) + g2(t)−G2, so that(4.36) is approximated bydxdt = b+H(x) + q(T ), (4.37)dbdt = −µ.H(x) = G1 +G2tanh(g3(x+ 1)) +G4x,q(T ) = g1(t)−G1 + g2(t)−G2,where Gj is the average value of gj(t) over one period for j = 1, 2, 4 and q(T )is an oscillatory function with mean zero and period of t = 1, with T = Ωtand Ω = 6.28. The computed trajectories shown in Figure 4.4 illustrate thevalidity of (4.37) as an approximation for (4.36).An analytical expression for the tipping point in the case q(T ) = 0 isfound as described in (4.4)-(4.5), based on a Taylor series of H(x) in (4.37)about x = 0. Thenbd = |H2|−1/3ad − bp for bp = H0 +H214|H2|, Hj = j!H(j)(0),is the tipping point for (4.37) for q(T ) = 0, and is compared below with thetipping point for q(T ) 6= 0. We note that H0 and H1 are small, but notidentically zero as in Section 4.1, since we are using an approximate system(4.37) rather than the full model from which the saddle-node bifurcationpoint (∆F0c, Ec) was determined.Next, we consider whether (4.37) corresponds to the high frequency forc-ing case as in Section 2.2.2, or the low frequency forcing case as in Section864.2. Energy Balance Sea Ice Model∆ F5 10 15 20E01020304050607080b-0.5 0 0.5 1 1.5x-2-1012345678Figure 4.4: Time series in both panels are shown for µ = .01. LEFT:Time series for E in (4.36) with oscillatory forcing terms gj (red dottedline) and without, that is, where the oscillatory terms gj are replaced bytheir averages Gj (black dash-dotted line). The critical value (∆F0c, Ec) ismarked with *. The solid blue curves correspond to dE/dt = 0 in (4.36)for the averaged model. RIGHT: Time series for x given in (4.37) withoscillatory forcing terms gj (red dotted line) and without, that is, where theoscillatory terms gj (black dash-dotted line), are replaced by their averagesGj . The solid blue curves correpond to b = −H(x) in (4.37). When rescaledusing the transformation (4.35), the oscillations for the x approximation areslightly smaller than those for the full E equation, although the differenceis negligible on this scale.874.2. Energy Balance Sea Ice Model2.3. For µ≪ 1 and Ω = 6.28, we would expect to be in the high frequencycase. However, since the amplitude of the oscillations is large, Section 2.2.3indicates that log(A)/ log(Ω) must be considered to determine if oscillationswith amplitude A correspond to high or low frequency forcing in a rescaledsystem. We find that 4/3 < log(max(q(T )))/ log(Ω) < 2 for (4.37) fallsin between the asymptotic ranges identified in Section 2.2.3 for the highand low frequency cases. Nevertheless, we use the concepts from the pre-vious sections to approximate the shift in the tipping point driven by theoscillations in (4.37).We begin with an outer approximation for x based on the approach fromSection 2.2.2 and Appendix B. Substituting a multiple scale approximation,x(µt, T ) = x0(µt, T ) + Ω−1x1(µt, T ) + . . . for T = Ωt into (4.37) we getx0T = 0,⇒ x0 = x0(t) x1T = b+H(x0) + q(T ). (4.38)The solvability condition (2.9) for x1 yields b+H(x0) = 0 so that the leadingorder contributions to x arex ∼ H−1+ (−b) + Ω−1Q(T ) , Q′(T ) = q(T ), (4.39)where H−1+ (−b) > 0 corresponds to the stable branch for x > 0 and Q(T )has zero average. The leading order approximation (4.39) is composed ofoscillations Ω−1Q(T ) about the branch of the averaged system b = −H(x).At this stage we appeal to the results shown for the canonical model inSection 2.2.2, noting that we must adapt that approach since the amplitude884.2. Energy Balance Sea Ice Modelof the oscillatory forcing is large. If the amplitude was O(1), we would writethe solution in the multiple scale form x = Ω−1Q(T )+Y (t, T ) analogous tothe local expansion in Section 2.2.2. This form is an approximation com-posed of oscillations near x = 0 with Y a small correction to be determined.However, since the oscillations are large for this application, we expect thatthe outer approximation breaks down away from x = 0. In fact, if we exam-ine the outer approximation (4.39), we see that it reaches values of x thatare on the unstable lower branch shown in Figure 4.4. We can approximatethe value of x where this occurs by identifying a value b∗ and correspondingx∗ = H−1+ (−b∗) on the stable branch of b = −H(x) satisfyingH−1+ (−b∗) + Ω−1 min(Q(T )) = H−1− (−b∗), (4.40)where H−1− (−b∗) is on the unstable branch shown in Figure 4.4. Then weadapt the approach from Section 2.2.2 by writing the solution in the multiplescale form that has oscillations about x∗ with Y a correction,x = x∗ +Ω−1Q(T ) + Y (T, t). (4.41)We could also rescale b, Y and the time scale t with µ. However, we obtainthe same result without this rescaling, so for simplicity we do not intro-duce such rescaled variables here. Substituting the multiple scale expression(4.41) into (4.37) yields894.2. Energy Balance Sea Ice ModelΩYT = −Yt + b+H(x∗ +Ω−1Q+ Y ) (4.42)∼ −Yt + b+H(x∗ +Ω−1Q) +H ′(x∗ +Ω−1Q)Y+ H′′(x∗ +Ω−1Q)2 Y2. (4.43)In contrast to the approximation used in (D.1) for the Morris-Lecarmodel or in (4.4), we have not replaced H(x) with a polynomial in x beforesubstituting (4.41). The reason that we avoid this step is that the forcingq(T ) is relatively large in this case, so we must consider H(x) over a largerrange of x values. Over this range, it is not possible to approximate H(x)accurately with a Taylor series about a single value of x. Rather we havesubstituted (4.41) directly into H(x) and then expanded about Y = 0,assuming that Y is a small correction. Then the solvability condition (2.9)for (4.43) yields the equation for YYt = b+(limL→∞1L∫ L0H(x∗ +Ω−1Q(T ))dT)+(limL→∞1L∫ L0H ′(x∗ +Ω−1Q(T ))dT)Y+(limL→∞1L∫ L0H ′′(x∗ +Ω−1Q(T ))2 dT)Y 2= b+H0 +H1Y +H2Y 2. (4.44)From (4.40) we get x∗ ≈ 2.2. Then we obtain the tipping point from the904.2. Energy Balance Sea Ice Modelsingularity for (4.44), using the results from (4.4) and getbtip = |H2|−1/3ad − bQ for bQ = H0 +H214|H2|. (4.45)In Figure 4.5 we compared this approximation with the numerical resultsin terms of the original variables E and ∆F0. As would be expected for asystem forced by oscillations with a large amplitude to frequency ratio, theadvance of the tipping point is dominated by the oscillations.ν0 0.01 0.02 0.03 0.04 0.05∆ F02468101214161820Figure 4.5: Comparison of the tipping points ∆F0 obtained from the asymp-totic approximation (4.45) (black dotted lines), compared with the numer-ically obtained values from the full model for E (4.36) (blue o’s) and fromthe approximate model for x (4.37) (red *’s) by the 2nd order Runge-Kuttamethod. The lower values correspond to the tipping point for the averagedsystem with gj(t) replaced with Gj , while the upper curve corresponds tothe case where the oscillations are not replaced with their averages.914.2. Energy Balance Sea Ice ModelThe fact that we are able to obtain an analytical expression for the tip-ping point allows us to explore the impact of parameters on the bifurcationstructure of the model. This information is valuable in quantitatively ob-taining the impact of parameter values on transitions between states. Morespecifically, we can explore the existence or prominence of hysteresis betweendifferent states as a function of the parameters of a periodically forced sys-tem. Within the context of a large scale atmosphere-ocean global climatemodel studied in [8], the question was considered whether there is hysteresisor not in the state transitions of (4.32) for E > 0 and E < 0. There itwas seen that the hysteresis predicted by lower dimensional models is lessprominent and in certain contexts not observed in the large scale compu-tations. The question of how the bifurcation structure could change withparameters was explored computationally in [15] for a simplified version of(4.32). For example, by recomputing the bifurcation diagram for differentparameter values, [15] shows a loss of the hysteresis loop with increased con-tributions to longwave radiation from FT (t), and thus a loss of hysteresisobserved for a varying bifurcation parameter. In the context of tipping, theloss of hysteresis corresponds to shifts in the tipping points so that regionsof bistability are no longer observed. The value of the analytical expressionfor the tipping point (4.45) is in determining shifts in the tipping locationas a function of a parameter, rather than having to recompute the entirebifurcation structure through repeated simulation of the full model over apotentially large range of parameter values.In Figure 4.6 we demonstrate how this loss of hysteresis is obtaineddirectly from the analytical expression (4.45) for the tipping point. The924.2. Energy Balance Sea Ice Model2 4 6 8 1051015202530c mlH ml∆ F 00 1 2 351015202530FT amplitude∆ F 0−2 −1 0 1 251015202530FT shift∆ F 00.3 0.4 0.5 0.6 0.7 0.851015202530hα∆ F 0Figure 4.6: Analytical approximation (blue circles) and numerical approxi-mation (red ∗’s) for the value ∆F0 corresponding to tipping as a function ofdifferent parameters: the ratio of E to T (t, E) in (4.34) (UPPER LEFT),a multiplicative factor for the amplitude of the oscillations FT (UPPERRIGHT), a shift in the mean of FT (LOWER LEFT), and the ice thicknessparameter describing transitions in the albedo (LOWER RIGHT). In theUPPER RIGHT figure, the approximation is not valid for tipping pointswith values ∆F0 > 31. This corresponds to larger amplitudes of FT , forwhich there is no overlap in the range of ∆F0 for upper and lower attractingbranches of E, indicating no hysteresis for varying ∆F0. Then determiningx∗ from (4.40) is no longer valid for the asymptotic approximation. In thatcase the numerical approximation asymptotes to ∆F0 ≈ 31.934.2. Energy Balance Sea Ice ModelUPPER RIGHT panel shows that increasing the amplitude of FT (t) causesearlier tipping. For tipping points increased to ∆F0 ≈ 31, there is a lossof hysteresis, since tipping occurs for approximately the same value of ∆F0for which there is a transition from the lower solution branch for E (notshown) to the upper branch in Figure 4.4. The approximation (4.45), basedon (4.40), indicates this behavior for the amplitude of FT increased above afactor of about 2.7, as also observed in [15]. For this increase in the oscil-lations there is no overlap in b for the upper and lower attracting branchesfor the solution x, and no intermediate unstable branch. Then determiningx∗ from (4.40) is no longer valid, signaling the loss of bistability and of hys-teresis. For the same reason, the numerical approximation asymptotes tothe value of ∆F0 ≈ 31, as shown in Figure 4.6(UPPER RIGHT). From thenormalized model (4.36) we can also identify other parameters like cmlHmland the mean of FT that can drive a considerable shift in the tipping point.The LEFT panels of Figure 4.6 show the effect of varying these parameters,which is not as large as for the increased amplitude of FT , as also observedin [15].The variation of the tipping point with hα is also shown. For this model,approximating the tipping location using (4.40)-(4.45) is based on the exis-tence of an upper solution branch. Then Figure 4.6(LOWER LEFT) showsthe approximation for values hα where x∗ in (4.40) can be determined. Theupper branch used in (4.40) no longer exists for hα sufficiently small, as alsoseen in [15].Summary of results: The method discussed in Chapter 2 can be adaptedfor the case with the amplitude of the periodic forcing larger than its fre-944.3. The Zero-dimensional Climate Modelquency, as was done for the energy balance sea ice model. To extract therelationship of parameters like the drifting rate and the frequency of the os-cillation, it is important to normalize the system (4.32) to a standard formwhere the dependent variable and bifurcation parameters take O(1) valuesand the saddle-node bifurcation point is located at the origin. Instead offocusing on the local region near x ∼ 0, we looked at the region near b = b∗where the leading order approximation of the normalized system reachesthe unstable equilibrium. Within the transformed model we can identifyparameter ranges where advanced tipping is likely to occur, by relating theresults of the canonical model to the rescaled parameters. The analyticalexpression of the location of the tipping point can be used to explore theinfluence of different parameters in the sea ice model, and further to explainthe change of the bifurcation structure [15].4.3 The Zero-dimensional Climate Model4.3.1 The modelThe zero-dimensional climate model is developed based on the globally av-eraged energy balance between the longwave emission and incoming solarradiation, including a linear albedo feedback [22]. The model is describedas follows [22, 29]:dT¯dt =1c{−ǫ0σT¯4 + 14 µ¯I0(1− a+ bT¯ )}, (4.46)where T¯ is the average surface temperature of the ocean on earth, ǫ0σT¯ 4 is954.3. The Zero-dimensional Climate Modelthe net outgoing longwave radiation described by a Stefan-Boltsmann law,14 µ¯I0(1 − a + bT¯ ) is the net incoming solar radiation, and c is the thermalinertia representing a thermal storage, estimated as 1.5×108 [kg/(K·s2)]. Forthe outgoing longwave emission, ǫ0 = 0.69 denotes the effective emissivity.The parameter σ is called a Stefan-Boltsmann constant and roughly equalsto 5.67×10−8 [J/(K4 ·m2 ·s)]. For the incoming shortwave radiation, a = 2.8and b = 0.009 [K−1] define a linear albedo feedback together. An upperbound for a− bT¯ is given as 0.75, to avoid complex or negative temperaturevalues. The parameter I0 = 1.36 × 103 [ (W/m2)] is the solar constant andµ¯ > 0 represents the variability of I0 [22].Theoretically, the system (4.46) has at most 4 equilibria because ofthe order of the polynomial on the right-hand side of (4.46). Consideringthe physical meaning of T¯ as the average surface temperature, we neglectthe complex roots. If we take µ¯ as the bifurcation parameter, there ex-ists a saddle-node bifurcation at the bifurcation point (T¯e = 4(a−1)3b , µ¯e =45ǫ0σ(a−1)333b4I0 ). For µ¯ > µ¯e, there exist a stable steady state (upper branch)and an unstable one(lower branch) [22], similar to the bifurcation diagramof the normalized system shown in Figure 4.7.If µ¯ is no longer a constant but slowly decreases to the critical value µ¯ewith a constant drifting rate 0 < v¯ ≪ 1, (4.46) becomesdT¯dt =1c{−ǫ0σT¯4 + 14 µ¯I0(1− a+ bT¯ )},dµ¯dt = −v¯,T¯ (0) = T¯0, µ¯(0) = µ¯0,(4.47)964.3. The Zero-dimensional Climate Modelwhere µ¯ = µ¯(t) means that it changes with time t, µ¯0 > µ¯e, and T¯0 is inthe neighbourhood of the stable branch of (4.46). The system (4.47) hasno equilibrium because µ¯ is slowly decreasing. The trajectory of (4.47) isclose to the stable steady state of (4.46), and decreases exponentially for µ¯is O(v¯2/3) smaller than µ¯e [32]. The abrupt transition from the stable stateto the unstable one, can be treated as a tipping point. The smaller that v¯is, the closer the tipping point is to the bifurcation point.If a noise term is introduced into (4.47), an early tipping is likely tohappen, depending on the amplitude of the noise [38, 41]. As a first step,we take an additive white noise. In order to consider the relative sizes of theparameters, we need to rescale (4.47) before the noise is introduced. Becausethe parameters take values from O(10−8) to O(108), terms like ǫ0σT¯ 4 areproducts of very small and very large quantities. We rescale the systemrelative to the bifurcation point and write in terms of quantities that arecomparable to O(1). We rescale time by τ = D · t and introduce two newparameters A = ǫ0σDc T¯ 3e and B =I0bDc . We normalize (4.47) by assumingT (τ) = 6A · T¯ (t)−T¯eT¯e and µ =3AB8 (µ¯− µ¯e). The system becomesdT =(c0T 4 + c1T 3 − T 2 + c2µT + µ)dτ + ǫdWτ ,dµdτ = −v,T (0) = T0, µ(0) =3AB8 (µ¯0 − µ¯e) = µini,(4.48)where c0 = − 163A2 , c1 = −19A , c2 = 23A , v = 3AB8D · v¯, ǫ denotes the amplitudeof the additive noise and Wτ represents the standard Wiener process. Tofocus on the early tipping before the bifurcation parameter reaches the crit-ical value, as shown in Figure 4.7, we consider µ ≥ 0 and ignore the possible974.3. The Zero-dimensional Climate Modelbifurcation structure for µ < 0.For ǫ = 0 and v = 0, (4.48) has a saddle-node bifurcation at (µe = 0, Te =0), which may have symmetric branches or not, depending on parametersc0, c1 and c2. Figure 4.7 shows the bifurcation structures of (4.48) withconstant µ (v = 0) and no noise (ǫ = 0) for zero and non-zero parametersc0, c1, and c2. For D = 12×107 [41], we have c0 = −0.4731, c1 = −1.1233,and c2 = 6.7396, which lead to an asymmetric bifurcation diagram and thestable branch of the bifurcation diagram is above that of the symmetriccase (c0 = c1 = c2 = 0), while the unstable branch is much flatter than thesymmetric case.4.3.2 Scaling relationship of the drifting rate and theamplitude of the noiseLet p(T, τ) be the PDF of T at time τ . The Fokker-Planck Equation (FPE)corresponding to (4.48) ispτ = −[(c0T 4 + c1T 3 − T 2 + c2µT + µ) · p]T +ǫ22 · pTT ,p(T, 0) = p0(T ).(4.49)Similar to Section 3.1, we introduce variables asv = ǫβ, µ = µini − vτ = ǫαs, T = ǫα/2ξ. (4.50)Then (4.49) becomes984.3. The Zero-dimensional Climate Model0 0.05 0.1 0.15 0.2−0.8−0.6−0.4−0.200.20.40.60.8µTFigure 4.7: Early tipping induced by relatively large noise (blue dotted)compared with the simulation of the slowly varying equilibrium solutionof the system (4.48) with ǫ = 0 (red dashed) and the bifurcation diagramof (4.48) with v = 0 (black solid and dash-dotted for stable and unstablebranches). Magenta solid and dash-dotted lines represent the stable and un-stable branches of the canonical model (3.1), which has the same bifurcationdiagram with (4.48) for c0 = c1 = c2 = 0. The amplitude of the noise isǫ = 0.05 and the drifting rate of the bifurcation parameter is v = 1.25×10−4.994.3. The Zero-dimensional Climate Model−ps = ǫ32α−βp[2ξ − ǫα/2(c2s+ 3c1ξ2)− ǫα4c0ξ3]+ ǫ 32α−βpξ[(ξ2 − s)− ǫα/2(c2sξ + c1ξ3)− ǫαc0ξ4] + ǫ2−βpξξ2 . (4.51)We assume that p(ξ, s) has the WKB form of the solution (3.5) [6].Substituting (3.5) into (4.51), we get−ψs = ǫ32α−β+γ [2ξ − ǫα/2(c2s+ 3c1ξ2)− ǫα4c0ξ3]+ ǫ 32α−βψξ[(ξ2 − s)− ǫα/2(c2sξ + c1ξ3)− ǫαc0ξ4]+ ǫ2−β−γψ2ξ2 + ǫ2−β ψξξ2 . (4.52)Similar to the analysis in Section 3.1, we substitute γ = 2 − 32α into(4.52) and get−ψs = ǫ32α−β[(ξ2 − s)− ǫα2 (c1ξ3 + c2sξ)− ǫαc0ξ4]ψξ+ ǫ2−β[2ξ − ǫα2 (3c1ξ2 + c2s)− ǫα4c0ξ3]+ ǫ 32α−β · 12(ψξ)2 + ǫ2−β · 12ψξξ. (4.53)Comparing (4.53) to (3.11), we notice that there are terms in O(ǫ2α−β),O(ǫ 52α−β), O(ǫ2−β+α2 ), and O(ǫ2−β+α) which add a complication to the fol-lowing analysis.1004.3. The Zero-dimensional Climate ModelBy substituting the general expansion (3.12) and (3.14) in (4.53) we havedifferent correction terms contribute to each order depending on the valueof α.For 0 < α < 45 , we haveO(ǫ 32α−β) :[(ξ2 − s)− ǫα2 (c1ξ3 + c2sξ)− ǫαc0ξ4]· ψ0ξ = −12(ψ0ξ)2,(4.54)O(ǫ 32α−β+δj ) :[(ξ2 − s)− ǫα2 (c1ξ3 + c2sξ)− ǫαc0ξ4]ψjξ + ψ0ξψjξ= Hj(ψ0, ψ1, · · · , ψj−1), j = 1, 2, · · · . (4.55)Unlike (3.7), −ǫα2 (c1ξ3+c2sξ)ψ0ξ−ǫαc0ξ4 in (4.54) contributes to the leadingorder O(ǫ 32α−β) because its order is lower than O(ǫ2−β). But similar to (3.7),(4.54) yields a quasi-stationary solution asψ0 = −2[(13ξ3 − sξ)− ǫα2 (14c1ξ4 + 12c2sξ2)− ǫα 15c0ξ5]+ f0(s), (4.56)where f0 is an undetermined function of time.To simplify the analysis, it is reasonable to ignore the term ǫ 52α−β(−c0ξ4ψξ)and ǫ2−β+α(−4c0ξ3) in the higher order case because T 4 is too small to makesignificant difference to the system (4.48) when T is close to 0, which is alsoconfirmed by the numerical simulations. By substituting (4.56) into theleft-hand side of (4.55), we have1014.3. The Zero-dimensional Climate Model[ǫα2 (c1ξ3 + c2sξ)− (ξ2 − s)]ψjξ = Hj(ψ0, ψ1, · · · , ψj−1), j = 1, 2, · · · .(4.57)Similar to the analysis of (3.16) and (3.17), the secularity condition of (4.57)is to eliminate the singularity at the stable equilibrium of (4.48).For j = 1, 2, we haveO(ǫ2−β) :[ǫα2 (c1ξ3 + c2sξ)(ξ2 − s)]ψ1ξ =− ǫβ−2ψ0s, (4.58)O(ǫ4−β− 32α) :[ǫα2 (c1ξ3 + c2sξ)− (ξ2 − s)]ψ2ξ =− ǫβ−2ψ1s −12(ψ1ξ)2− 12ψ1ξξ. (4.59)Substituting (4.56) into (4.58) and change variables back to T and µ by(4.50), we getψ1T =ǫβ−2(2T + c2T 2 + ǫ32αf0µ)T 2 − µ− c1T 3 − c2µT= ǫβ−2(2T + c2T 2 + ǫ32αf0µ)−c1(T − T1)(T + T2)(T + T3), (4.60)where T1, −T2, −T3 are the roots of the equation T 2−µ− c1T 3− c2µT = 0.The minus sign in front of T2 and T3 indicates there are two negative rootsdepending on the value of c1 and c2. In order to find f0, we eliminate thesingularity at T = T1 which is the stable equilibrium of (4.48) with v = 0and c0 = 0. We have1024.3. The Zero-dimensional Climate Modelf0µ = ǫ−32α[−(2T1 + c2T 21 )]⇒ f0(µ) = ǫ−32α∫ µ0−2T1(z)− c2T1(z)2dz.(4.61)By substituting (4.61) into (4.60), we haveψ1 = ǫβ−21c1(−T3 + T2){[2 + c2(T1 − T2)] ln(T + T2)− [2 + c2(T1 − T3)ln(T + T3)]}+ f1(µ), (4.62)where f1(µ) is determined in the following steps.Substituting (4.56), (4.61) and (4.62) into (4.59) and applying the secu-larity condition similar to (3.17), we getf1(µ) =∫ µ1[ǫ2−β(−12ψ11TT −12ψ211T )− ψ11µ]|T=T1(z)dz, (4.63)where ψ11 denotes ψ1(T, µ) without f1(µ).Therefore, (4.56), (4.61), (4.62) and (4.63) provide the asymptotic ap-proximation of PDF p(T, µ), with a normalization coefficient C0:p(T, µ) ∼ C0 · exp{ 1ǫ2−3α/2ψ0(T, µ) + ψ1(T, µ)}, (4.64)where C0 is estimated byC0 =[∫ Tp(z, µ(τ))dz]−1. (4.65)1034.3. The Zero-dimensional Climate ModelFor α ≥ 45 and β + 32α < 4, we can apply similar procedures like (4.53)-(4.65) to get asymptotic results. However, it would be balance less terms inthe leading order and higher orders since −ǫα2 (c1ξ3 + c2sξ) in (4.54) is nolonger a correction of O(ǫ 32α−β). Hence, to simplify the analysis, we use thesame asymptotic approximation as (4.64).4.3.3 PDF tail approximation for trajectories crossing theunstable branchSimilar to Section 3.3, after rescaling T = ǫη, we havepτ +µ0ǫ pη −12pηη = ǫ[η2pη + 2ηp −µ0ǫ c2(p+ ηpη) + ǫβ−2τpη]+ ǫ2[c1(3η2p+ η3pη + ǫβ−2c2τ(p+ ηpη))], (4.66)p(η, 0) =pqs(η, µ0), η < Ts/ǫwhere µ0 ∼ O(ǫ4/3), pqs is the asymptotic approximation (4.64).Substituting the expansion of p(η, τ) as (3.26), we find that the leadingorder contributions are the same as in (3.27) and the solution is the sameas (3.28). The higher order distributions are as follows:1044.3. The Zero-dimensional Climate ModelO(ǫ) : p1τ +µ0ǫ p1η −12p1ηη = f1(p0, η, τ)= η2p(0)η + 2ηp0 −µ0ǫ c2(p0 + ηp0η)+ ǫβ−2τp0η, (4.67)O(ǫ2) : p2τ +µ0ǫ p2η −12p2ηη = f2(p0, p1, η, τ)= η2p(1)η + 2ηp1 −µ0ǫ c2(p1 + ηp1η)+ ǫβ−2τp1η + c1(3η2p0 + η3p0η)+ ǫβ−2c2τ(p0 + ηp0η), (4.68)· · ·which have the solutions in the original variables as followspj(T, τ) =∫ τ0∫ ∞−∞1ǫ√2π(τ − s)exp[−(T − µ0τ − z)22ǫ2(τ − s)]· fj(p0, ..., pj−1(z, s), z, s)dzds, (4.69)where fj, j = 1, 2 are shown in (4.67) and (4.68).Similar to the analysis in Section 3.3, we introduce Td = −T2+2ǫ, belowwhich the approximation (4.69) breaks down, and Tm = −T2 + ǫ where wepatch (4.69) with (4.64).1054.3. The Zero-dimensional Climate Model4.3.4 Numerical results and comparison to the analyticalapproximationWe compare the PDF given by (4.64) and (4.69) with the results basedon numerical simulations by the Euler-Maruyama method in order to seethe possibility of the early tipping under different choices of parameters.Our approximation is used to estimate the probability of early tippingas P (T¯ < T¯unstable) =∫ T¯unstable−∞ p(T¯ , µ¯(t))dT¯ = 1 −∫∞T¯unstable p(T¯ , µ¯(t))dT¯ ,where T¯unstable represents the unstable equilibrium at µ¯ = µ¯(t), and differ-ent indicators of early tipping, for example, the skewness and variance ofthe PDF. By changing the values of c0, c1 and c2, we change the structureof the bifurcation diagram and explore whether it changes the probabilityof the early tipping.Figure 4.8 shows the PDF of (4.48) given by (4.64) for µ¯ is relatively farfrom the critical value µ¯e = 0.9655. The probability of the early tipping isnegligible. Under the same relationship (β = 3) and at the same value of thebifurcation parameter µ¯ = 1.7878, the PDF with larger noise (ǫ = 0.1) hasa heavier left tail, shown in Figure 4.8(RIGHT), which indicates a largerchance of early tipping while µ¯ keeps decreasing. For smaller noise (ǫ =0.05), even though the drifting rate v is 10 times smaller (corresponding tov¯ = 1.03 × 10−10 in (LEFT) and v¯ = 8.22 × 10−10 in (RIGHT)), the PDFis dense around its center, shown as Figure 4.8(LEFT), which indicates asmaller probability of early tipping. Figure 4.8 also shows the influence of thebifurcation structure. The system with c0 = c1 = c2 = 0 has a symmetricbifurcation structure, which is the same as the one for the system (3.1),1064.3. The Zero-dimensional Climate Model250 300 350 40000.0050.010.0150.020.0250.030.0350.04T¯p(T¯,µ¯)200 250 300 350 40000.0050.010.0150.02T¯p(T¯,µ¯)Figure 4.8: Comparison between different PDF approximations for µ¯ =1.7878 relatively far from its critical value µ¯e = 0.9655. LEFT: ǫ = 0.05,β = 3, v = ǫβ = 1.25 × 10−4. RIGHT: ǫ = 0.1, β = 3, v = ǫβ = 0.001. Bluesolid line represents the asymptotic approximation for system (4.48) givenby the WKB approximation (4.64). Black dashed curve is the asymptoticapproximation of the PDF for the symmetric case (c0 = c1 = c2 = 0) givenby the WKB approximation (3.23). Magenta dash-dotted curve is based onthe Gaussian type approximation (3.10). Red circle represents the numericalhistogram under 2.5 × 105 realizations.1074.3. The Zero-dimensional Climate Modelshown as the black dashed curves in Figure 4.8. At the same µ¯ and for thesame ǫ and β, the PDF of the system with the flatter stable branch (blackdashed) shows larger probability for the system to be near the unstableequilibrium T¯ < T¯e = 266.67, which indicates a larger probability of earlytipping for smaller µ¯ than the one with a steeper stable branch (blue solid).200 250 300 35000.0050.010.0150.020.0250.03T¯p(T¯,µ¯)µ¯= 1.2122150 200 250 300 35000.0050.010.0150.02T¯µ¯= 1.0478Figure 4.9: PDF at different µ¯ closer to its critical value µ¯e = 0.9655 wherea tipping is more likely to happen for the case ǫ = 0.05 and β = 2.5.Blue solid line represents the asymptotic approximation for system (4.48)given by the WKB approximation (4.64) and the tail approximation (4.69).Black dashed curve is the asymptotic approximation of the PDF for thesymmetric case (c0 = c1 = c2 = 0) given by the asymptotic approximation(3.23) and (3.30). Magenta dash-dotted curve is based on the Gaussiantype approximation (3.10). Red circle represents the numerical histogramunder 2.5×105 realizations. Blue dotted vertical line represent the unstableequilibrium for (4.48) with c0 = −0.4731, c1 = −1.1233, and c2 = 6.7396.Blue and black stars in Right represent T¯d and T¯s for the approximation ofthe left tail of the PDF and patching process by normalization.For µ¯ smaller (µ¯ = 1.21), the left tail of the PDF starts to cross the un-stable equilibrium and the approximation (4.64) cannot provide approxima-tion for T¯ below the unstable equilibrium T¯unstable. We combine the WKB1084.3. The Zero-dimensional Climate Modelapproximation (4.64) with the left tail approximation (4.69) by patching.The results are shown in Figure 4.9(LEFT). For ǫ = 0.05 and β = 2.5, theapproximation given by (4.64) and (4.69) agrees with the numerical simula-tions while the left tail of the PDF starts to cross the unstable equilibrium(T¯unstable = 243.98 with c0 = −0.4731, c1 = −1.1233, and c2 = 6.7396,shown as the blue dotted vertical line), indicating that the early tippingis likely to happen. For µ¯ even smaller (µ¯ = 1.05), the PDF shown inFigure 4.9(RIGHT) indicates a larger probability for T¯ below the unstableequilibrium T¯unstable = 251.01. Even though (4.69) cannot be used for T¯far below T¯unstable, the approximation of PDF given by (4.64) and (4.69)agrees with the numerical simulations for T¯ > T¯d, which indicates thatthe probability of early tipping at certain values of µ¯ can be estimated byP (T¯ < T¯unstable) = 1 − P (T¯ > T¯unstable) = 1 −∫∞T¯unstable p(T¯ , µ¯(t))dT¯ . Theapproximation of the probability of early tipping P (T¯ < T¯unstable) doesn’tmean the system (4.48) escapes from the stable state as soon as T¯ belowthe unstable equilibrium T¯unstable, but provides an over-estimation. If thesystem (4.48) is in the neighbourhood of the unstable equilibrium at µ¯, itescapes at smaller µ¯ due to the convection-diffusion process.Figure 4.10(LEFT) shows the probability of T¯ above T¯unstable as 1 −P (T¯ < T¯unstable) = P (T¯ > T¯unstable) =∫∞T¯unstable p(T¯ , µ¯(t))dT¯ at any µ¯ > µ¯efor ǫ = 0.05 and β = 2.5, for the symmetric canonical model (c0 = c1 =c2 = 0) and temperature model (4.48). The vertical dotted line shows thelocation where the quasi-stationary approximation (4.64) is no longer validand where we introduce the left tail approximation (4.69). For decreasing µ¯,P (T¯ > T¯unstable) decreases below 1, which indicates an increasing probabil-1094.3. The Zero-dimensional Climate Modelity that T¯ crosses the unstable equilibrium of the system at µ¯ and decreasesexponentially before µ¯e = 0.9655. The Gaussian approximation cannot cap-ture the value of µ¯ where the early tipping starts. The approximation givenby (3.23) and (3.30) based on the symmetric bifurcation diagram providesa much larger probability of early tipping.1 1.1 1.2 1.3 1.40.40.50.60.70.80.911.1µ¯P(T¯>T¯unstable)1 1.1 1.2 1.3 1.40.40.50.60.70.80.911.1µ¯P(T¯>T¯unstable)Figure 4.10: LEFT: P (T¯ > T¯unstable) estimated by (4.64) and (4.69) (Bluesolid) and estimated by numerical simulations (Black dashed with star mark-ers) for (4.48) with c0 = −0.4731, c1 = −1.1233, and c2 = 6.7396, estimatedby (3.23)(Red dash-dotted), numerical simulation (Black dashed) and Gaus-sian type approximation (Magenta dotted with circle markers) for (4.48)with c0 = c1 = c2 = 0. Vertical blue dotted line shows the region wherethe WKB approximation is no longer valid for all the values of T¯ (α > 43 ).Right: Comparison between the probability P (T¯ > T¯unstable) based on thesystem (4.48) with c0 = −0.4731, c1 = −1.1233, and c2 = 6.7396 (blue solidfor asymptotic approximation and black dashed with markers for numericalsimulations) for ǫ = 0.05 and v = 0.0005, and P (T¯ > T¯unstable) based onthe system (4.48) with c0 = 0, c1 = 1.1233, and c2 = −6.7396, estimatedby (4.64) and (4.69) (red dash-dotted) and numerical simulations (blackdashed).To see the influence of the asymmetry of the bifurcation structure on theprobability of early tipping, we compare the probability P (T¯ > T¯unstable)1104.3. The Zero-dimensional Climate Modelbased on system (4.48) with c0 = −0.4731, c1 = −1.1233, and c2 = 6.7396,and the system (4.48) with c0 = 0, c1 = 1.1233, and c2 = −6.7396 for ǫ =0.05 and v = 0.0005, shown in Figure 4.10(RIGHT). For c0 = 0, c1 = 1.1233,and c2 = −6.7396, the related bifurcation diagram has a flatter stable branchand steeper unstable one than the ones for c0 = −0.4731, c1 = −1.1233, andc2 = 6.7396. Under the same size of noise and drifting rate, the probabilityof early tipping P (T¯ < T¯unstable) at the same µ¯ becomes much larger if theshape of the stable equilibrium branch is flatter, and P (T¯ > T¯unstable) forthe case c0 = 0, c1 = 1.1233, and c2 = −6.7396 is smaller than the one forc0 = −0.4731, c1 = −1.1233, and c2 = 6.7396.Figure 4.11 shows the change of the probability for T¯ > T¯unstable if wedouble the noise (UPPER) or reduce the drifting rate to half (LOWER)from µ = 0.05 and v = 0.0005 in the form of probability (LEFT) and theconcentration region (RIGHT) bounded by the dashed lines in which theprobability of T is 95%. For ǫ = 0.1, the early tipping is likely to happenwhen µ¯ < 1.8, which is much earlier than the one for ǫ = 0.05, shown inFigure 4.11(UPPER LEFT). From the concentration region shown in Figure4.11(UPPER RIGHT), we notice that the larger noise induces a heavier lefttail of the PDF which crosses the unstable branch of the bifurcation diagrammuch earlier than the case with the smaller noise. For v = 0.00025, the earlytipping is likely to happen slightly earlier than the one for v = 0.0005, shownin Figure 4.11(LOWER LEFT). The concentration region shown in Figure4.11(LOWER RIGHT), indicates that the PDFs of the system (4.48) areidentical at each µ¯ until the left tail of the PDF with the slower drifting ratev = 0.00025 crosses the unstable branch of the bifurcation diagram slightly1114.3. The Zero-dimensional Climate Modelearlier than the one with the faster drifting rate v = 0.0005.Changing skewness of the distribution of time series data is claimed tobe a reliable warning signal of tipping [42, 43], and so is increasing vari-ance [43]. In [19], the skewness and variance of the stationary density werecomputed for a saddle-node system with an additive Gaussian noise and aconstant bifurcation parameter. The results based on the empirical station-ary density, which is the approximation of the stationary density functionbased on samples, show an increasing trend of the skewness while the bifur-cation parameter gets closer to the bifurcation point. The approximation ofthe skewness obtained from time series data doesn’t show an obvious trendof change compared with the skewness based on the empirical stationarydensity, but the sign of the skewness estimated from the time series doesshow the existence of the nonlinear terms in the system. Those studies arebased on the system with a constant bifurcation parameter. For a systemwith slowly varying bifurcation parameter like (3.1) and (4.48), we can ex-plore whether the change of the skewness and variance can be used as awarning signal by the PDF approximation (4.64) and (4.69). We estimatethe skewness and variance of the PDF for every µ¯ byTˆ =∫T¯ p(T¯ , µ¯(t))dT¯ , (4.70)σ2 =∫(T¯ − Tˆ )2p(T¯ , µ¯(t))dT¯ , (4.71)γ =∫(T¯ − Tˆ )3p(T¯ , µ¯(t))dT¯σ3 . (4.72)1124.3. The Zero-dimensional Climate Model1 1.2 1.4 1.6 1.80.20.30.40.50.60.70.80.911.1µ¯P(T¯>T¯unstable)1 1.2 1.4 1.6 1.8150200250300350400µ¯T¯1 1.2 1.4 1.60.30.40.50.60.70.80.911.1µ¯P(T¯>T¯unstable)1 1.2 1.4 1.6 1.8150200250300350400µ¯T¯Figure 4.11: LEFT: P (T¯ > T¯unstable) with different parameter values. Theblue solid (black dashed with stars) curve is for the case ǫ = 0.05, v =0.0005 given by (4.64) and (4.69) (numerical simulations). UPPER LEFT:The red dash-dotted (black dashed) line is for the case ǫ = 0.1, v = 0.005given by (4.64) and (4.69) (numerical simulations). LOWER LEFT: The reddash-dotted (black dashed) line is for the case ǫ = 0.05, v = 0.00025 givenby (4.64) and (4.69) (numerical simulations). RIGHT: The concentrationregion near the stable equilibrium bounded by T¯h and T¯l at different µ¯,which are defined by P (T¯ > T¯h) =∫∞T¯h p(T¯ , µ¯(t))dT¯ = 0.025, P (T¯ < T¯l) =∫ T¯h−∞ p(T¯ , µ¯(t))dT¯ = 0.025. The black solid (dash-dotted) line represents thestable (unstable) branch of the bifurcation diagram of (4.46). Blue dashedlines are for the case ǫ = 0.05, v = 0.0005. UPPER RIGHT: Red dashedlines are for the case ǫ = 0.1, v = 0.0005. LOWER RIGHT: Red dashedlines are for the case ǫ = 0.05, v = 0.00025.1134.3. The Zero-dimensional Climate Model1.2 1.3 1.4 1.5 1.6−0.8−0.7−0.6−0.5−0.4−0.3−0.2−0.10µ¯Skewnessγ (I)(II)(III)1.2 1.3 1.4 1.5 1.60100200300400500600µ¯Varianceσ2Figure 4.12: ǫ = 0.05, v = 0.0005. The skewness γ (4.72) and variance σ2(4.71) based on asymptotic approximation of the PDF. The Roman numeralscorrespond to cases with different values of c0, c1 and c2: Case (I) c0 =−0.4731, c1 = −1.1233, and c2 = 6.7396, estimated by (4.64) and (4.69);Case (II) c0 = c1 = c2 = 0, estimated by (3.23) and (3.30); Case (III) c0 = 0,c1 = 1.1233, and c2 = −6.7396, estimated by (4.64) and (4.69).1144.3. The Zero-dimensional Climate ModelFigure 4.12 shows the skewness and variance estimated based on p(T¯ , µ¯(t))given by (4.64) and (4.69) for different bifurcation structures. The minussign of γ represents the bias of the left tail of the PDF. We notice thatall three cases have the same qualitative behaviour of the skewness andthe variance. For both symmetric and asymmetric case, the magnitudeof the skewness based on the asymptotic approximation (4.64) and (4.69)increases when the bifurcation parameter is relatively large but slowly ap-proaches the critical value. And the magnitude of the skewness of the PDFstarts to decreases while µ¯ is closer to its critical value and continues to de-crease. Those behaviours of the skewness γ are qualitatively consistent withthe results based on the empirical stationary density in [19]. For all threecases, the variance keeps increasing while µ¯ decreases. Comparing Figure4.10(RIGHT) with Figure 4.12, we notice that for all three cases, the magni-tude of the skewness of the PDF increases before the system (4.48) reachesthe value of µ¯ where P (T¯ < T¯unstable) is no longer negligible. This increasingtrend of the skewness as a warning signal for the tipping is consistent withthe results in [42, 43].In [19], the approximation of the skewness obtained from time series datadoes not show an obvious trend of change compared with the skewness basedon the empirical stationary density. The values of the skewness from timeseries are much smaller than those from the empirical stationary density.And in [19], it says that only the sign of the skewness indicates the nonlinearterm in the system which induces tippings in the future. However, if wecompare the range of values of the skewness for the symmetric case c0 =c1 = c2 = 0 (Case (II) in Figure 4.12) to the results from time series in [19],1154.3. The Zero-dimensional Climate Modelthey are close to each other. The consistency between our results and theresults from time series in [19] shows that even though the change of theskewness is small due to the slowly varying bifurcation parameter, relative tothe skewness based on the empirical distribution in [19], it is not negligibleand it reveals the probability of early tipping.In Figure 4.12, we notice that among the three cases, the magnitude ofthe skewness γ for Case (I) is the largest at different µ¯ and the magnitudeof the skewness for Case (II) reaches its maximum at the largest µ¯. Basedon the results in [19, 42, 43], we would expect the earliest tipping in eitherCase (I) or Case (II). However, the probability P (T¯ > T¯unstable) in Figure4.10 shows the system of Case (III) has the largest probability of tippingat every µ¯, which is inconsistent with the possible conclusion that mightbe reached based on using only the skewness in Figure 4.12. The reasonwhy γ cannot reveal the probability of early tipping for different bifurcationstructures is that it is affected by many factors, like the mean of the PDFand the variance. Based on the WKB approximation of the PDF (4.64), weknow that the system (4.48) is most likely to stay in the neighbourhood ofthe stable equilibrium corresponding to the local minimum of the potentialwell, which is close to Tˆ in (4.70). If the system (4.48) has a steeper stablebranch of the bifurcation diagram, at the same µ¯, the PDF has a smallervariance since the majority of the PDF is close to Tˆ , which leads to a largervalue of the skewness γ since the variance σ2 is shown as the denominatorin (4.72). Therefore, the system of Case (I) has the largest magnitude of theskewness and the smallest variance in Figure 4.12, but the system of Case(I) has the smallest probability of the early tipping at every bifurcation1164.3. The Zero-dimensional Climate Modelparameter µ¯. To consider the influence of the asymmetry of the bifurcationstructure on the early tipping, it is important to use the approximation ofthe PDF (4.64) and (4.69), instead of the tipping indicators like skewnessand variance.Generally speaking, the asymmetry of the bifurcation structure is in-duced by the coefficients c0, c1 and c2 in (4.48), of which the most impor-tant factor is c2 as the coefficient of term µT . If c2 = 0, the system can berenormalized around the bifurcation point and the new bifurcation diagramis symmetric. For c2 > 0, the stable branch of the bifurcation diagram issteeper and the unsable branch is flatter than those for c2 = 0, which in-dicates a smaller probability of early tipping. For c2 < 0 which providesa flatter stable branch compared to the symmetric case c2 = 0, the sys-tem (4.48) has a larger chance of early tipping. The importance of c2 alsoshows that instead of directly comparing the historical data with canonicalmodel (3.1), it is necessary to explore the relationship between the bifur-cation parameter and the system variable, in other words, the existence ofthe coefficient c2. The asymptotic approximation of the PDF (4.64) and(4.69) can provide useful information about the asymmetry of the bifurca-tion diagram based on the historical data, and how early the tipping couldhappen.Summary of results: The WKB method and multiple-scale analysis wedeveloped in Chapter 3 can be applied on systems similar to the canonicalmodel (3.1) which have a saddle-node bifurcation, a slowly varying bifur-cation parameter and an additive noise, for example, the zero-dimensionalclimate model with an asymmetric bifurcation diagram. The comparison be-1174.3. The Zero-dimensional Climate Modeltween the Gaussian-type approximation and the WKB approximation showsthat the asymmetry of the PDF reveals the possibility of the early tipping.The comparison of P (T¯ > T¯unstable) for the system with different shapesof the bifurcation diagrams shows that the asymmetry of the bifurcationdiagram plays an important role on how early the system can escape fromthe stable state. For example, the system with a flatter stable branch has alarger probability of escaping while the systems reach the same value of thebifurcation parameter. In [42, 43], the skewness of PDF from the canoni-cal system (3.1) with constant bifurcation parameter has been used as anindicator of tipping. In [19], a comparison of skewness between the the-oretical results from canonical model with constant bifurcation parameterand the ones from a time series has been discussed, but not for the systemwith slowly varying bifurcation parameter. The influence of the drifting rateof the bifurcation parameter is included in the skewness estimated by theasymptotic approximation of the PDF (4.64) and (4.69), and so is the asym-metry of the bifurcation structure. For systems with different shapes of thebifurcation diagram, the indicators of early tipping, such as the skewnessor the variance of the PDF show the same qualitative behaviour, indicatingthe existence of tipping. But there’s not enough quantitative informationto tell how early the tipping could happen for the systems with differentbifurcation structures under the same drifting rate and random forcing. Itis important to compare the PDF given by the time dependent data and theasymptotic approximation to understand the bifurcation structure of thesystem and how early the tipping can happen.118Chapter 5Conclusion and Future WorkWe have analyzed the factors that determine the location of a tipping pointin the canonical model for a saddle-node bifurcation with bifurcation param-eter a that varies at a slow rate µ ≪ 1 and forced by a periodic oscillationwith amplitude A and frequency Ω. We see that the main contributions thatdetermine the location of the tipping point are different for low and highfrequencies. In order to get analytical expressions to identify the key factorsin the location of the tipping point, we apply the methods of multiple scalesand matched asymptotic expansions to obtain analytical expressions for thecombined effect of Ω, A, and µ. The cases of low and high frequency in theoscillatory forcing naturally require the consideration of different relation-ships between the time scales of the oscillation and the slow drifting rate.The size of the amplitude also plays an important role.For a high-frequency forcing, there are two competing contributions tothe location of the tipping point relative to the static bifurcation point ata = 0. The first contribution is a delay due to the slow drifting rate µ ofthe bifurcation parameter a, well-known in a variety of studies without pe-riodic forcing [31, 32]. The second contribution is an advance proportionalto the square of the amplitude to frequency ratio of the oscillation, ob-119Chapter 5. Conclusion and Future Workserved also in the case where the bifurcation parameter is constant (µ = 0).These contributions reflect that the trajectory in this case is approximatedby the linear combination of the slowly varying equilibrium solution foundfor A = 0 in [32] and an additive contribution that takes the form of theoscillatory forcing with a scale factor. Near the tipping point, the trajec-tory can be expressed in terms of an Airy function on the slow time scale,with the oscillations averaged on this time scale shifting the singularity thatcorresponds to the tipping point. The resulting expressions for these con-tributions then give an expression for the location of the tipping point thatvaries continuously with the parameters. Whether the combined effect re-sults in an advance or delay compared to a = 0 depends on the magnitudeof µ2/3 relative to A2/Ω2. In particular, we see noticeable advances of thetipping point for combinations of larger amplitudes with frequencies thatare not very large.Advances in the tipping point due to periodic forcing are even moreapparent for lower frequencies. For Ω ≪ 1, the trajectory is dominated bya nonlinear combination of the slowly drifting bifurcation parameter andthe oscillatory forcing, since they give contributions that are of the sameorder of magnitude. Then the location of the tipping point depends notonly the amplitude, frequency, and drifting rate, but also on the phase ofthe oscillation. Tipping occurs near values of a where the trajectory isnear x ∼ 0, where a local analysis reveals that the change of concavityfor x ∼ 0 provides an analytical expression for the tipping point. Thisexpression captures the discontinuous dependence of the tipping locationon parameters such as the amplitude. These discontinuities follow directly120Chapter 5. Conclusion and Future Workfrom the fact that tipping occurs only for values of a where the trajectory isnear zero: a slight change in parameters can change the approach to zero ofthe trajectory, postponing tipping to a later point at which the trajectorycrosses zero. This shift in the tipping point is on the order of one period ofthe low frequency oscillation.We have analyzed the factors that determine the probability of the es-cape in the canonical model for a saddle-node bifurcation with bifurcationparameter a that varies at a slow rate µ ≪ 1 and affected by an additivewhite noise with an amplitude ǫ. In order to find the PDF of the systemby solving the related FPE, we use the WKB method for a relatively farfrom the critical value 0. The leading order WKB approximation is relatedto the potential well of the system as a quasi-stationary approximation.The approximation including higher order terms shows that the PDF ei-ther approaches to 0 or it becomes a singularity at the unstable equilibriumwhere the approximation is no longer valid. For a close to 0, to capture thenon-zero PDF near the unstable equilibrium, we introduce a multiple-scaleanalysis, considering the system near the stable and unstable equilibriumseparately, by the WKB approximation and a convection-diffusion process.With the approximation of the PDF at any value of the bifurcation param-eter, we are able to explore the relationship between the early tipping andthe amplitude of the noise, the drifting rate of the bifurcation parameter,even the asymmetric structure of the saddle-node bifurcation diagram.For a far from the critical value 0, the WKB approximation can beapplied for β > 2, which reflects the relationship between the drifting rate µof the bifurcation parameter and the amplitude of the noise ǫ. When β < 2,121Chapter 5. Conclusion and Future Workwhich is the same as µ > ǫ2, the Gaussian-type approximation given by[34] is applicable, which indicates the sample paths of the system is likelyto follow the slowly varying equilibrium solution of the system and no earlytipping happens. When β > 2, which is the same as µ < ǫ2, the leadingorder WKB approximation is driven by the potential well of the system atevery value of the bifurcation parameter. For smaller a, the potential wellbecomes flatter, and the PDF of the system has a heavier left tail whichindicates the system has larger probability to escape from the stable stateif a is closer to the critical value 0. The behaviour of the system separatedby the relationship µ = ǫ2 is consistent with the conclusion of [28] based ona sample-path analysis.For smaller a where the stable and unstable equilibrium are both closerto 0 and no longer far from each other, the WKB approximation isn’t ap-plicable near the unstable equilibrium. Therefore we introduce a local anal-ysis of the system near x = 0 and a = 0, which turns out to be a seriesof convection-diffusion equations. By finding the solution asymptoticallyand patching it with the WKB approximation near the stable equilibriumthrough normalization, we find the approximation of the PDF of the sys-tem which provides the probability of the early tipping for any a > 0. Thecombined PDF shows that for a approaching 0, non-zero PDF rises near theunstable equilibrium and the probability of the system below the unstableequilibrium is no longer negligible, which indicates that it is possible for thesystem to escape from the stable state for those a > 0.The analysis on the Morris-Lecar model with a slowly increasing inputcurrent and an additive periodic external input shows that our method de-122Chapter 5. Conclusion and Future Workveloped in Chapter 2 is applicable to the system with a bifurcation structuresimilar to a saddle-node bifurcation, a slowly varying bifurcation parameterand an additive periodic forcing, but in higher dimensions. For the highfrequency oscillation, the location of the tipping point is determined by thelocal bifurcation structure near the bifurcation point, the slowly drifting bi-furcation parameter, the amplitude and the frequency of the periodic forcing.For the low frequency case, tipping occurs near values of b where the trajec-tory is near x ∼ 0, which is determined by not only the varying bifurcationparameter, but also the amplitude, the frequency and the initial conditionof the external oscillatory forcing. The discontinuity of the location of thetipping point related to the amplitude of the periodic forcing is captured bythe concavity analysis similar to the one in Chapter 2. The analysis on theenergy balance sea ice model shows that the method discussed in Chapter2 can be adapted for the case with the amplitude of the periodic forcinglarger than its frequency. Instead of focusing on the local region x ∼ 0, welooked the region near b = b∗ where the leading order approximation of thenormalized system reaches the unstable equilibrium. In both applications,the influence of the oscillations, and thus the choice of the appropriate ap-proximation, depend on certain scaling relationships between the driftingrate of the bifurcation parameter, amplitude, and frequency. These areidentified within the context of the canonical model where the dependentvariable and bifurcation parameters take O(1) values and the saddle-nodebifurcation point at the origin. Thus an important step is normalization ofthe model to this standard form, from which the parameter relationshipscan be extracted. Within the transformed model one can already identify123Chapter 5. Conclusion and Future Workparameter ranges where advanced tipping is likely to occur, by relating theresults of the canonical model to the rescaled parameters.We applied the WKB method and multiple-scale analysis on the zero-dimensional climate model which has a saddle-node bifurcation, a slowlyvarying bifurcation parameter and an additive noise to analyze the effect ofthe factors like the amplitude of the noise, the drifting rate and the struc-ture of the bifurcation diagram on the probability of early tipping. Thecomparison between the Gaussian-type approximation and the WKB ap-proximation shows that the asymmetry of the PDF reflects the possibilityof the early tipping. The comparison of P (T¯ > T¯unstable) for the systemwith different shapes of the bifurcation diagrams shows that the asymmetryof the bifurcation structure plays an important role on how close the bifur-cation parameter can be to the bifurcation point for the system to escapefrom the stable state. At the same value of the bifurcation parameter, thesystem with a flatter stable branch has a larger probability of escaping. Forsystems with different shapes of the bifurcation diagram, the indicators ofearly tipping, such as the skewness or the variance of the PDF show thesame qualitative behaviour, but there is not enough quantitative informa-tion to tell which system has larger probability to escape under the samedrifting rate and random forcing. It is essential to identify the relationshipbetween the bifurcation parameter and the system variable, in other words,the existence of the factors which affect the shape of the bifurcation diagram.There are a number of remaining questions about the influence of oscil-latory and noisy forcing on a system with a slowly drifting bifurcation. Hereand in [19, 28] the case of additive forcing has been considered, but in gen-124Chapter 5. Conclusion and Future Workeral the oscillatory and noisy forcing can be multiplicative or parametric, ascan be seen explicitly in the model considered in [16] without averaging overseasonal periodicity. In addition, the combination of noisy and oscillatoryforcing should be of interest, as they appear together in applications. Forexample, for low frequency oscillations, we see that the location of the tip-ping point depends on the phase, and we would expect that the noise causesphase shifts. The question then remains when and if stochastic phase shifts,in combination with periodic forcing, play a significant role in the tippinglocation. A related area is the phenomenon of rate induced tipping, whereexceeding a critical drifting rate can result in early tipping [4]. The effectof noisy or oscillatory forcing on these types of transitions has not beenanalyzed.The shift of the tipping point observed for different types of oscillationsindicates that we must consider their effect on transitions via a saddle-nodebifurcation. For example, a common underlying structure for a hysteresisloop consists of two branches for stable states with a bistability region andloss of stability via saddle-node bifurcations. If oscillations are included thatadvance the tipping points, the size of the hysteresis loop can be reduced,or the hysteresis can be removed completely.More generally in a slow-fast system where the slow variable plays therole of a bifurcation parameter for the rest of the system, any oscillatoryforcing can change the time to transition when a saddle-node bifurcation isembedded in the systems structure. Complex phenomena such as burstingdynamics or mixed-mode oscillations typically appear in slow-fast systemsand involve combinations of transitions via different types of bifurcation125Chapter 5. Conclusion and Future Workwith a slow variable viewed as a bifurcation parameter. For phenomena ofthis type in which a saddle-node bifurcation is present, the results we havecan quantify the influence of oscillations on the length of the intervals wheresteady states are attracting. In addition to the saddle-node bifurcationconsidered here, there are other types of bifurcations that are potentiallysensitive to combinations of slowly drifting parameters and oscillatory andnoisy forcings. This sensitivity has been studied in other contexts, includ-ing different types of bifurcations [7, 28, 34], although not necessarily fromthe point of view of considering competitive or combined effects of bothoscillations and noise on transitions or bifurcations.126Bibliography[1] Maykut G. A. and N. Untersteiner. Some results from a time-dependentthermodynamic model of sea ice. J. Geophys Res, 76:1550–1575, 1971.[2] Sutera A. On stochastic perturbation and long-term climate behaviour.Quarterly Journal of the Royal Meteorological Society, 107(451):137–151, 1981.[3] D.S. Abbot, Silber M., and Pierrehumbert R.T. Bifurcations leadingto summer arctic sea ice loss. J. Geophysical Res., 116:D19120, 2011.[4] R. Vitolo Ashwin P., S. Wieczorek and P. Cox. Tipping points in opensystems: bifurcation, noise-induced and rate-dependent examples in theclimate system. Phil. Trans. R. Soc. A, 370:1166–1184, 2012.[5] T. Erneux Baer S.M. and J. Rinzel. The slow passage through a hopfbifurcation: Delay, memory effects, and resonance. SIAM J. AppliedMath., 49:55–71, 1989.[6] C.M. Bender and S.A. Orszag. Advanced Mathematical Methods forScientists and Engineers. McGraw-Hill, New York, 1978.[7] P. Borowksi and R. Kuske. Characterizing noisy mixed mode oscilla-tions in neuronal models. Chaos, 20:043117, 2010.127Bibliography[8] Armour K. C., I. Eisenman, E. Blanchard-Wrigglesworth, K.E. Mc-Cusker, and C.M. Bitz. The reversibility of sea ice loss in a state-of-the-art climate model. Geophysical Research Letters, 38:L16705, 2011.[9] Meisel C. and C. Kuehn. Scaling effects and spatio-temporal multileveldynamics in epileptic seizures. PLos ONE, 7(2):e30371, 2012.[10] Morris C. and H. Lecar. Voltage oscillations in the barnacle giant mus-cle. Biophys J., 35:193–213, 1981.[11] He D. and D. D. Earn. Epidemiological effects of seasonal oscillationsin birth rates. Theoretical Polulation Biology, 72:274–291, 2007.[12] T. Erneux and J.-P. Laplante. Jump transition due to a time-dependentbifurcation parameter in the bistable ioadate-arsenous acid reaction. J.Chem. Phys., 90:6129–6134, 1989.[13] Strogatz S. H. Nonlinear Dynamics and Chaos. Westview Press, Read-ing, MA., 1994.[14] A. Hohl, H.J.C. van der Linden, G. Goldstein R. Roy, F. Broner, andS.H. Strogatz. Scaling laws for dynamical hysteresis in a multidimen-sional laser system. Phys. Rev. Lett., 74:2220–2223, 1995.[15] Eisenman I. Factors controlling the bifurcation structure of sea iceretreat. Journal of Geophysical Research, 117:D01111, 2012.[16] Eisenman I. and J. S. Wettlaufer. Nonlinear threshold behavior duringthe loss of arctic sea ice. Proc. Natl. Acad. Sci. U. S. A., 106(1):28–32,2009.128Bibliography[17] Grasman J. and O. A. van Herwaarden. Asymptotic Methods for theFokker-Planck Equation and the Exit Problem in Applications. SpringerNew York, New York, 1999.[18] Merryfield W. J., M. M. Holland, and A. H. Monahan. Multiple equi-libria and abrupt transition in arctic summer sea ice extent. ArcticSea Ice Decline: Observations, projections, mechanisms, and implica-tions, Bitz, C. M., Deweaver, E. T., Tremblay, L. B., Eds. AmericanGeophysical Union:151–174, 2008.[19] Sieber J. and J. M. Nonlinear softening as a predictive precursor to cli-mate tipping. Philosophical Transactions of the Royal Society A: Math-ematical, Physical and Engineering Sciences, vol. 370, issue 1962:1205–1227, 2012.[20] Zhu J., R. Kuske, and T. Erneux. Tipping points near a delayed saddlenode bifurcation with periodic forcing. SIAM J on Applied DynamicalSystems, Submitted.[21] P. Jung, G. Gray, and R. Roy. Scaling law for dynamical hysteresis.Phys. Rev. Lett., 65:1873–1876, 1990.[22] Fraedrich K. Structural and stochastic analysis of a zero-dimensionalclimate system. Quart. J. R. Met. Soc., 104:461–474, 1978.[23] J. Kevorkian and J.D. Cole. Multiple Scale and Singular PerturbationMethods, volume 114 of Appl. Math. Sciences. Springer New York, NewYork, 1996.129Bibliography[24] Dai L., D. Vorselen, K. S. Korolev, and J. Gore. Generic indicators forloss of resilience before a tipping point leading to population collapse.Science, 336:1145, 2012.[25] Lenton, H. Held T. M., E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf,and H. J. Schellnhuber. Tipping elements in the earth’s climate system.Proc. Nat. Acad. Sci. USA, 105(6):1786–1793, 2008.[26] Pal M., A. K. Pal, S. Ghosh, and I. Bose. Early signatures of regimeshifts in gene expression dynamics. Phys. Biol, 10:036010, 2013.[27] Scheffer M, Carpenter S. R., Lenton T.M., Bascopmpte J., Brock W.,Dakos V., van de Koppel J., van de Leemput I. A., Levin S. A., van NesE. H., Pascual M., and Vandermeer J. Anticipating critical transitions.Science, 338:344–348, 2012.[28] Berglund N., in C. Laing Barbara Gentz, and G. J. Lord (Edts).Stochastic methods in neuroscience. Oxford University press, pages65–93, 2010.[29] US National Institute of Standards and Technology. Codata value:Stefan-boltzmann constant. The NIST Reference on Constants, Units,and Uncertainty, 2011.[30] Dafilis M. P., F. Frascoli, J. McVernon, J. M. Heffernan, and J. M. Mc-Caw. The dynamical consequences of seasonal forcing, immune boostingand demographic change in a model of disease transmission. Journal ofTheoretical Biology, 361:124-32:doi: 10.1016/j.jtbi.2014.07.028., 2014.130Bibliography[31] Mandel P. and T. Erneux. The slow passage through a steady bifurca-toin: delay and memory effects. J. Stat. Phys., 48:1059–1070, 1987.[32] Haberman R. Slowly varying jump and transition phenomena associ-ated with algebraic bifurcation. SIAM J. Applied Math., 37(1):69–106,1979.[33] Haberman R. Elementary Applied Partial Differential Equations withFourier Series and Boundary Value Problems. Prentice Hall, UpperSaddle River, 1998.[34] Kuske R. Probability densities for noisy delay bifurcations. Journal ofStatistical Physics, Vol. 96, Issue 3-4:797–816, 1999.[35] Herrmann S. and P. Imkeller. The exit problem for diffusions with timeperiodic drift and stochastic resonance. Ann. Appl. Probab. 15, No.1A:39–68, 2005.[36] Lee S. Study on the onset bifurcations of a morris-lecar neuron undera periodic current. Journal of the Korean Physical Society, 50:346–350,2007.[37] Lee S. Bifurcation analysis of a mode-locking structure in a stronglyforced morris-lecar neuron. Journal of the Korean Physical Society,52:11–16, 2008.[38] Thompson J. M. T. and J. Sieber. Climate tipping as a noisy bifur-cation: a predictive technique. IMA Journal of Applied Mathematics,76:27–46, 2011.131[39] K.K. Tung and J. Zhou. Using data to attribute episodes of warm-ing and cooling in instrumental records. Proc. Nat. Acad. Sci. USA,110(6):2058–2063, 2013.[40] Anagnostopoulou V. and T. Ja¨ger. Non-autonomous saddle-node bi-furcations: random and deterministic forcing. J. Diff., 253(2):379–399,2012.[41] Dakos V., M. Scheffer, E. H. van Nes, V. Brovkin, V. Petoukhov, andH. Held. Slowing down as an early warning signal for abrupt climatechange. Proc. Nat. Acad. Sci. USA, 105(38):14308–14312, 2008.[42] Guttal V. and C. Jayaprakash. Changing skewness: an early warningsignal of regime shifts in ecosystems. Ecol Lett. 2008 May, 11(5):450–460, 2008.[43] Guttal V. and C. Jayaprakash. Spatial variance and spatial skewness:leading indicators of regime shifts in spatial ecological system. TheorEcol, pages DOI 10.1007/s12080–008–0033–1, 2008.[44] Wikipedia. Morris-lecar model — Wikipedia, the free encyclope-dia, 2015. Retrieved from http://en.wikipedia.org/wiki/Morris-Lecarmodel, on 14-April-2015.[45] Y. Li Yu N. and R. Kuske. A computational study of spike time re-liability in two types of threshold dynamics. Journal of MathematicalNeuroscience, 3:11, 2013.132Appendix AThe Higher Order Analysisfor µ = 0 and Ω ≫ 1We continue the asymptotic analysis in Section 2.2.1 after (2.11) for O(Ω−2)and O(Ω−3), yieldingO(Ω−2) :x2T + x1t = −2x0x1,⇒ x2T = R2(T, t) = −v1t − 2x0v1 + 2Ax0cos(T ) (A.1)O(Ω−3) :x3T + x2t = −x21 − 2x0x2,⇒ x3T = R3(T, t) = −x2t − x21 − 2x0x2, (A.2)By applying the solvability condition (2.9) to (A.1) and (A.2) and solvingthe resulting equations, we haveO(Ω−2) : v1t = −2x0v1, (A.3)x2T = 2Ax0cos(T ), ⇒ x2 = 2√aAsin(T ) + v2(t), (A.4)O(Ω−3) : v2t = −A22 − 2√av2, (A.5)133Appendix A. Local Analysis for High Frequency ForcingThe attracting solutions of (A.3) and (A.5) are their stable equilibria, whichare v1 = 0 and v2 = − A24√a . After substituting these equilibria, (2.11) and(A.4) into (2.6), we get (2.12) in Section 2.2.1.In order to find the local behaviour of (2.1) with µ = 0, we substitutea = Ω−2b into (2.5) and getxT +Ω−1xt = Ω−3b− Ω−1x2 +Ω−1Asin(T ). (A.6)The contributions at each order areO(Ω0) : x0T = 0, ⇒ x0 = x0(t), (A.7)O(Ω−1) : x1T + x0t = −x20 +Asin(T ),⇒ x1T = −x0t − x20 +Asin(T ), (A.8)O(Ω−2) : x2T + x1t = 0, ⇒ x2T = −x1t, (A.9)O(Ω−3) : x3T + x2t = b− x21, ⇒ x3T = −x2t + b− x21.(A.10)Substituting (A.7) into (A.8), applying the solvability condition (2.9),and solving the resulting equations, we havex0t = −x20, ⇒ x0 = 0, (A.11)x1T = Asin(T ), ⇒ x1 = −Acos(T ) + v1(t). (A.12)Substituting (A.12) into (A.9) and applying the solvability condition (2.9),134Appendix A. Local Analysis for High Frequency Forcingwe havev1t = 0, ⇒ v1 = d, ⇒ x2T = 0, ⇒ x2 = x2(t), (A.13)where d is constant and needs to be determined. Substituting (A.12) and(A.13) into (A.10) and applying the solvability condition (2.9), we havex2t = b−A22 − d2, ⇒ d2 = b− A22 , (A.14)which provides the equilibrium solution of x2. If b − A22 < 0, x2t is alwaysnegative and there is no attracting solution.By substituting (A.11), (A.12), (A.13) and (A.14) into (2.6), we havethe leading order approximation (2.13) in Section 2.2.1.135Appendix BThe AsymptoticApproximation of theSolution for 0 < µ≪ 1 andΩ ≫ 1 in the Outer RegionWe substitute (2.15) into (2.14), yieldingO(1) : x0T = 0, ⇒ x0 = x0(τ),O(µλ) : x1T = a− x20 +Asin(T ). (B.1)Applying the solvability condition (2.9) of (B.1) and solving the resultingequations yieldsx20 = a, ⇒ x0 =√a,x1T = Asin(T ), ⇒ x1 = −Acos(T ) + v1(τ).136Appendix B. Local Analysis for Low Frequency ForcingThe higher order contributions of (2.14) depend on the value of λ. The twocases are:For 0 < λ ≤ 1, O(µ2λ) : µ1−λx0τ + x2T = −2x0x1. (B.2)For λ > 1, O(µ1+λ) : x0τ + x2T = µλ−1 · (−2x0x1). (B.3)In either case, we substitute x0 =√a in (B.2) or (B.3) and apply thesolvability condition (2.9), to getx0τ = µλ−1(−2x0v1) ⇒ v1 = µ1−λ ·14a.Then the asymptotic approximation of the attracting solution is the result(2.16) of Section 2.2.2.137Appendix CThe Local Approximation forΩ ≪ 1, near amWriting (2.26) in terms of the local variables η, ξ, and Am by (2.33), wehave−µξη = am[1− cos(cµη)] + µη −1c sin(cµη) − µA1amAmcos(cµη)− µ A1cAmsin(cµη)− µ2ξ2. (C.1)Substituting the expansion ξ ∼ ξ0 + µξ1 + µ2ξ2 + · · · , into (C.1), yieldsO(µ) : −ξ0η = −A1amAm, ⇒ ξ0 =A1amAmη + c0, (C.2)O(µ2) : −ξ1η =12amc2η2 − A1Amη − ξ20 ,⇒ ξ1 = −(12c2am −A21a2mA2m)η33 +A1Am(1 + 2c0am)η22 + c20η + c1,(C.3)which provides the local approximation (2.36). We find c0 by writing (2.36)in terms of x and an intermediate scaled variable a − am = Kµq, for q < 1138Appendix C. Local Analysis for Amplitude Thresholdand K an O(1) positive constant, to getx = A1amAm· Kµq + µc0 − (12c2am −A21a2mA2m)K3µ3q−13+ A12Am(1 + 2c0am)K2µ2q + µ1+qc2K + µ2c1 + · · · , (C.4)and match (C.4) with the Taylor expansion of the outer approximation (2.27)around am written in terms of Kx = µq ·√c2am2 K+ µ1−q · 12K (1−A1amAmc√2am) +O(µ2q, µ2−3q). (C.5)We take 13 < q < 12 so that the outer approximation (C.5) is valid. Bymatching the lower order terms in (C.4) and (C.5), we find c0 in the formof c0 = − 12am + C1 + C2. This form is convenient for comparing withthe ξ-intercept −1/(2am) of the concavity threshold (2.35), reducing thiscomparison to determining the sign of C1 + C2, given byC1 = µ−qAmKamA1(√c2am2 −A1amAm) + µ1−3q Am2K3amA1(1− A1amAmc√2am),C2 = µq−1Am3amA1(12c2am −A21a2mA2m).We find c0 for A > Am only, since Ac > Am by the definition of Am and am.If A1amAm <c2Am2A1 and A > Am, both C1 and C2 are positive and c0 > −12am .Then the local trajectory (2.36) is in the concave up region as assumed forAc > A > Am and a ≤ am (η ≤ 0). If A1amAm >c2Am2A1 , both C1 and C2 arenegative. In that case the local trajectory (2.36) is in the concave down139Appendix C. Local Analysis for Amplitude Thresholdregion for a ≤ am (η ≤ 0), resulting in tipping near a = am.140Appendix DThe Local Approximation forw in the ML Model withΩ ≫ 1The equations for the corrections z1 and W1 arez1T ∼ − µ1/3z0s + b+ (h00 + h01(X + z0) + h02(X + z0)2)− gK(X + z0 +D)W0, (D.1)W1T = −µ1/3W0s +(κ(0) + κ′(0)(X + z0) +κ′′(0)2 (X + z0)2)·(w∞(0) + w′∞(0)(X + z0) +w′′∞(0)2 (X + z0)2 −W0), (D.2)141Appendix D. Local Analysis for ML Modelwhere we have kept linear and quadratic terms in X and z0, analogous to(2.20). We apply the solvability condition as in (4.12), yielding for W0W0 = ( limL→∞L−1∫ L0(κ(0) + κ′(0)(X + z0) + κ′′(0)(X + z0)2/2)·(w∞(0) + w′∞(0)(X + z0) + w′′∞(0)(X + z0)2/2) dT )/(limL→∞L−1∫ L0(κ(0) + κ′(0)(X + z0) + κ′′(0)(X + z0)2/2) dT)≈ w∞(0) +W00 +W01z0 +W02z20 (D.3)W00 =w′∞(0)κ′(0) + w′′∞(0)κ(0)/2W0dA22Ω2W01 =w′∞(0)κ(0)W0d− κ′(0)W00W0dW02 = W00(2Ω2A2 −κ′′(0)2W0d+ 2κ′(0)2W 20d)− w∞(0)κ′(0)κ(0)W 20dW0d = κ(0) +A2κ′′(0)4Ω2where we have used a Taylor series about z0 to obtain a second orderpolynomial in z0 as an approximation for W0. Substituting (D.3) in theequation for z0 written in terms of the original time scale yields (4.16).142Appendix EGaussian Type Solution forthe Canonical ModelFor α > 23β and β < 2, we change variables byǫη = x− F (t), p(x, t)dx = pˆ(η, t)dη, (E.1)where F (t) is the solution of the deterministic system (ǫ = 0) in order tolook for a solution of (3.1) of the formpˆ(η, t) =√g(t)exp−g(t)2[η − 2√gn(t)]2, (E.2)where gn(t) is related to the mean and g(t) is related to the change ofvariance over time.Substituting (E.1) into (3.2), we havepˆt = ǫ−1pˆη{F ′(t)− [µ−F 2(t)]}+12 pˆηη +2pˆF (t)+ 2ηpˆηF (t)+ ǫ(2ηpˆ+ η2pˆη).(E.3)The terms contribute to O(ǫ−1) vanish because of the definition of F (t).143Appendix E. Gaussian ApproximationThe leading order of (E.3) is given by O(ǫ0),pˆt ∼12 pˆηη + 2F (t)(pˆ + ηpˆη). (E.4)By substituting (E.2) into (E.4), we have g(t) satisfying− g′2 =12g2 − 2gF (t), (E.5)which givesg(t) = e4Q2 ∫ t0 e4Q(z)dz+cg≡ gn2 ∫ t gn(z)dz ,Q(t) =∫ t F (z)dz,(E.6)and (3.10) in Section 3.1.144
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivs 2.5 Canada
Graduate
Asymptotic behaviour of tipping points in non-autonomous systems with random and periodic forcing
Text