Noise-induced tipping under periodic forcing: Preferred tipping phase in a non-adiabatic forcing regime

In this paper, we focus on a simple stochastic differential equation for which we can readily control the characteristic time scales of the problem. Specifically, we consider

Here, X_t ∈ ℝ is a stochastic process parameterized by time t ≥ 0, W_t is a standard Wiener process, and σ > 0 denotes the noise strength. We choose the parameters ɛ, α, A > 0 so that the deterministic problem has two stable periodic solutions separated by an unstable one [represented by solid and dashed curves in Fig. 1(a)]. The parameter ɛ represents a ratio of the characteristic relaxation time of the flow to the period of the forcing and arises when we nondimensionalize time by the forcing period. Throughout the paper, the “lower” and “upper” stable periodic solutions of $x˙=f(x,t)/ε$ are denoted by $xl∗$ and $xu∗$⁠, respectively, and the “middle” unstable one is denoted by $xm∗$⁠, where $xl∗(t)<xm∗(t)<xu∗(t)$ for t ∈ ℝ. The presence of noise introduces a mechanism for transition from one stable periodic solution to the other. Figure 1(a) shows some example realizations of (1). Because we are interested in the phase of tipping events relative to the forcing, rather than the absolute timing of expected tipping, we often wrap noisy realizations of (1) to a cylinder; see Fig. 1(b).

In Fig. 1(a), we present a schematic diagram of different ɛ and σ regimes in which noise-induced transitions differ qualitatively. Comparing insets in Fig. 1(a), we see that the value of ɛ in (1) can strongly impact the geometry and separation of the periodic solutions. This timescale ratio, in turn, may affect when and how the sample paths cross the unstable periodic solution. The adiabatic regime (regime II) in Fig. 1(a), where ɛ ≪ 1, has been well studied.^8,10,12–17 In this case, the relaxation time to stable periodic orbits is fast compared to the forcing period. Consequently, noise-induced transitions are “instantaneous” and occur with high probability when the separation between the stable and unstable periodic orbits is minimal and when an associated potential barrier height is lowest.

In contrast, regime I of Fig. 1(a), obtained with $ε=O(1)$ and otherwise the same parameters, has been the focus of only a small number of studies.^18–20 In this regime, the transition occurs more slowly, especially after crossing the unstable periodic solution, and the separation between the deterministic solutions does not vary significantly over the forcing period. The distribution of tipping events may display more than one narrow peak when crossing the unstable limit cycle.¹⁹ Moreover, in this regime, where noise is not too large, the impacts of noise and deterministic drift may act on comparable time scales. These observations suggest challenges with identifying a preferred tipping phase outside of the adiabatic regime. Indeed, while formal matched asymptotic expansions of solutions to the Fokker–Planck equation yield quantitative information about the distribution of tipping events when crossing the unstable orbit,¹⁹ little is known about what deterministic quantities govern when the most probable paths depart from the lower stable orbit. Such deterministic indicators are important for determining when the system is most susceptible to noise-induced transitions.

In the noise-dominated regime [regime III in Fig. 1(a)], randomness overpowers the deterministic dynamics and tipping samples are extremely noisy, with deterministic features largely washed out. In the fast forcing regime, marked as regime IV in Fig. 1(a), where the timescale ratio ɛ is large, the relaxation time of the dynamics is slow compared to the forcing period and transitions can take a long time compared to the forcing period, possibly making the notion of the preferred tipping phase moot.

As shown in Fig. 2(a), we let t_l, t_m, and t_u denote the times when a sample path departs from $xl∗$⁠, first crosses $xm∗$⁠, and arrives at $xu∗$⁠, respectively. They are defined as

Because we are interested in the phase of tipping, relative to the forcing, rather than the absolute time, we will typically mod these times by the forcing period T and speak of the “tipping phases” that then derive from t_l, t_m, and t_u. In Figs. 2(b) and 2(c), we present the distributions of t_l, t_m, and t_u modded by T = 1. We take deterministic parameters in regime I and consider noise intensity σ = 0.15 for Fig. 2(b) and σ = 0.3 for Fig. 2(c). We notice that for both values of σ, the histograms of t_l have a more prominent peak than the histograms of t_m and t_u. Moreover, the location of the peak in t_l does not shift with the change in σ the way that the peak in t_m does. Indeed, we find that the histograms for values of σ ∈ [0.15, 0.4] (not shown here) give consistent results; the peak of t_l broadens but does not shift significantly with increasing noise strength in this range. We will identify t_l as a robust “tipping phase” in this intermediate regime I. In addition to looking for conditions on ɛ and σ under which we can identify a robust “tipping phase,” we aim to investigate what features of the deterministic dynamics select such a phase when it exists.

In this paper, we study noise-induced transitions between stable periodic solutions in the intermediate forcing regime [regime I shown in Fig. 1(a)]. We use the framework of “exit problems,”²¹ where the most probable path of escape, from one domain of attraction to another, is typically determined by taking the limit σ → 0. In this setting, “most probable paths” (also called “optimal paths”) are often approximated by calculating the critical points of either the Freidlin–Wentzell (FW)^{18,19,21–23} or the Onsager–Machlup (OM) rate functional, which include an additional Ito correction term to the FW functional.^24,25 In the autonomous setting, the use of the FW functional can be rigorously justified using the theory of large deviations,^21,26 and for autonomous drift with a single time scale, the FW and OM functionals yield essentially identical most probable paths.^27,28 Efficient computational techniques for computing most probable paths range from direct minimization algorithms,^29,30 ordered upwind methods applied to the corresponding Hamilton–Jacobi equations,³¹ finite difference schemes that exploit the Hamiltonian structure of the Euler–Lagrange equations,³² and numerical continuation methods.³³ For nonautonomous systems with periodic forcing, the analysis is nuanced since the additional time scale inherent to the periodic forcing may introduce discrepancies between minimizers of the FW and OM functionals. Moreover, for our specific problem, consideration of other small parameters in the problem, such as the timescale ratio ɛ, may be important in resolving discrepancies.

The rest of the paper is organized as follows. In Sec. II, we formulate the problem using a path integral approach. We describe the FW and OM functionals that we consider, and how the parameters ɛ and σ enter the variational problem for computing optimal paths. Section III contains the main contributions of this paper. Focusing on regime I where ɛ = O(1), we perform numerical simulations over a range of parameters in order to identify a robust definition of the preferred tipping phase. We find that it is associated with the time of initial departure from the lower stable periodic solution, as suggested by the histograms of Figs. 2(b) and 2(c). It occurs at a phase of the forcing when the flow near this stable periodic solution changes, momentarily, from negative to positive, i.e., from downwards to upwards in Fig. 2(a). This phase is determined by an intersection of the periodic solution with the nullclines f(x, t) = 0 of the deterministic system. We conclude in Sec. IV with a discussion of our results, limitations of the method, and some avenues for future research. We further numerically investigate the most probable transition phase in a conceptual Arctic sea ice model due to Eisenman and Wettlaufer,⁷ which has a strong seasonal forcing and possesses bistable parameter regimes. We use insights from our case study to conjecture that this model is in regime IV of Fig. 1(a), the fast forcing case.

To simplify the notation, throughout the paper, we use $x˙$ to represent time derivative dx/dt and subscript to denote partial derivative, e.g., V_x(x, t) : = ∂V/∂x.

We first discuss the role of the parameters in the deterministic problem

where V(x, t) is a potential associated with the deterministic dynamics f(x, t)/ɛ, which is the drift term in (1)

The parameter A determines the forcing strength and α breaks a reflection symmetry of the potential. The timescale ratio ɛ controls how (un)stable the periodic solutions $x∗$ of (2) are; their Floquet multipliers, $ρ(x∗)$⁠, depend exponentially on 1/ɛ

Figure 3 summarizes the dependence of the deterministic dynamics of (2) on the parameters. The phase diagram in Fig. 3(b) shows that (2) admits four qualitatively distinct parameter regions for α > 0; these are characterized by different configurations of the nullclines f(x, t) = 0, in the (t, x)-phase plane. In this paper, we focus on regions R₁ − R₃ of the (A, α)-parameter plane where, for a given ɛ, there are two stable periodic solutions. (There is only one stable periodic solution in R₄.) The location of the saddle-node bifurcations of Fig. 3(a) depends on ɛ as indicated. We describe these three parameter regions and their nullclines as follows:

• The nullclines in R₁, defined by $0<α+A<2/(33)$⁠, consist of three curves that partition extended phase space into four regions with vertical flow directions that alternate in sign. A standard Poincaré analysis proves existence of two stable periodic solutions for all ɛ values in R₁.

• Region R₂, defined by $2/(33)<α+A$ and $|α−A|<2/(33)$⁠, has two nullcline curves that divide the phase space into three regions.

• In R₃, defined by $0<α<A−2/3$⁠, the phase space is partitioned by the nullcline into two regions.

Figure 3 highlights another ɛ-dependent feature of the dynamics. When ɛ ≪ 1 in (2), the periodic solutions are well-approximated by the nullclines f(x, t) = 0, as seen in Fig. 3(c) obtained for ɛ = 0.01. In contrast, for the moderate case ɛ = 0.5 of Fig. 3(d), there are significant deviations between the periodic solutions and the nullclines, leading to a pronounced intertwining of the (blue and black) curves in that phase portrait. This indicates intervals of time when the flow is upwards, say, along the lower stable solution, i.e., toward the upper unstable solution, and so on. These intervals become more pronounced in R₂ and R₃, where, provided ɛ is large enough, three periodic solutions may exist even though there are fewer than three curves comprising the nullcline; in these regions, the double well structure of V(x, t) is temporarily lost during certain phases of the forcing.

We use the path integral formulation to identify a possible preferred phase for tipping. For t_f > t₀ ≥ 0, let $P={x∈C([t0,tf];R):x(t0)=x0\ and\ x(tf)=xf}$ denote the set of paths connecting a point (t₀, x₀) on $xl∗$ and a point (t_f, x_f) on $xu∗$⁠. If $x∈P$ is differentiable, then the probability of the tipping trajectories lying in an infinitesimally small tubular neighborhood of $x,Pσ[x]$ satisfies

where d_W [x] is the Wiener measure,^25,34 and the proportionality reflects the missing normalization prefactor. Heuristically, this formulation can be derived by considering an equally spaced partition t₀ = t₁ < t₂ < ⋅ ⋅ ⋅ < t_n = t_f of width Δt, calculating the probability that a realization passes through a sequence of intervals [a_i, b_i] at time t_i, and taking the limit |b_i − a_i| → 0 and Δt → 0. Derivations following this approach can be found in the physics literature,^25,34 and we include our version of the formal derivation in  Appendix A.

The probability of a tipping event lying in a small tubular neighborhood of a tipping path $x∈P$ is then determined by integrating $Pσ[x]$ over this neighborhood. Similar to Laplace’s method for the asymptotic expansion of exponential integrals, we expect as σ → 0 that the dominant contribution to this integral is concentrated about local minimizers of the argument of the exponential in the definition of (5). With this as motivation, we define the “most probable path” or “optimal path” connecting points (t₀, x₀) and (t_f, x_f) to be minimizers of the Onsager–Machlup (OM) functional $Iσ:A↦R$ given by

which can be expressed as the sum of the Freidlin–Wentzell (FW) rate functional I_FW²¹ and an additional integral I_OM2. Here, the admissible set is $A:={x∈H1([t0,tf];R):x(t0)=x0,x(tf)=xf}$ and $L:([t0,tf],R,R)↦R$ denotes the Lagrangian. We note that for our model problem, the integrand of I_OM2 is simply f_x = 1 − 3x², so that this is typically a positive contribution to $Iσ$ when the paths are near the unstable solution $xm∗$ and a negative contribution when the paths are near the stable solutions $xl∗$ and $xu∗$⁠.

Since $L(t,x,x˙)$ is convex in $x˙$⁠, to prove the existence of a minimum in $A$⁠, it is sufficient to show that I_σ is coercive with respect to the H¹ norm. The next theorem establishes this fact (as $∥x˙∥L2→∞,Iσ→∞$⁠); the existence of a minimizer is then a corollary following standard techniques from the direct method of the calculus of variations.³⁵ 

Local minimizers of (6) satisfy the Euler–Lagrange equations:

These equations constitute a second order nonlinear boundary value problem in x on the domain [t₀, t_f] and can be rewritten in the Hamiltonian form, for time-varying Hamiltonian function

as

Here, the “momentum” $Ψ=x˙−ε−1f(x,t)$ measures the deviation from the deterministic flow. The FW functional in (6) can then be expressed in terms of Ψ as $IFW[x(t)]=∈tt0tfΨ(x,t)2dt$⁠.

The path integral formulation assumes that the noise strength σ is small. If σ is sufficiently small compared to the timescale ratio ɛ, such that σ²/ɛ ≪ 1, the naive regular perturbation in σ of (9), at leading order in σ, is

In this system, Ψ = 0 forms an invariant submanifold foliated by solutions of the deterministic dynamics $x˙=ε−1f(x,t)$⁠. Thus, for σ → 0, optimal paths are well-approximated by the heteroclinic orbits of (9) connecting deterministic solutions.^36,37 Moreover, (10) is the Hamiltonian form of the Euler–Lagrange equations for the FW rate functional I_FW. Since I_OM2 is independent of $x˙$⁠, minimizers of I_σ will converge uniformly as σ → 0 to a minimizer of I_FW; see  Appendix B.

The remainder of the paper addresses the intermediate regime I in Fig. 1(a), where ɛ is $O(1)$⁠. In this case, there is no small parameter to exploit to simplify the analysis. We investigate the contribution I_OM2 to the functional (6) in determining optimal paths, which are then validated by Monte Carlo simulations. This section concludes with a parameter study of (1), based on the path integral method using the Onsager–Machlup functional (6).

To solve the second order boundary value problem (7), we introduce a gradient flow $us=−δIσ[u]δu$ associated with $Iσ$⁠. Specifically, introducing the artificial “time” s ≥ 0, we consider solutions $u(s,t):R+×[t0,tf]↦R$ to the following partial differential equation:

where x_g(t) is a sufficiently smooth initial guess of the most probable path. To numerically solve (11), we use MATLAB’s built-in parabolic initial-boundary value solver pdepe.³⁸ 

Stationary solutions of (11) correspond to solutions of the Euler–Lagrange equations (7). If we let u(s, t) denote a solution to (11), we find [after integrating by parts and applying (11)] that

where we have suppressed the arguments in f and u in order to make the notation compact. Hence, I_σ[u(s, ⋅ )] is a monotonically decreasing function in s. Furthermore, by Theorem 1, I_σ[u(s, ⋅ )] is bounded from below, and, thus, the pointwise limit lim _s→∞I_σ[u(s, ⋅ )] exists and $lims→∞ddsIσ[u(s,⋅)]=0$⁠. I_σ is a Lyapunov function on $A$⁠, and, by Corollary 2, lim _s→∞u(s, t) = x(t), where x(t) is a local minimizer of the OM functional and, hence, a solution to the Euler–Lagrange equations.³⁹ Convergence to a steady state is assessed by the value of $δIδu$⁠, i.e., when the Euler–Lagrange equations are approximately satisfied.

For the gradient flow, we impose a starting position at $x0=xl∗(t0)$ and a final position at $xf=xu∗(tf)$⁠. Figure 4 presents results for six forcing periods between t₀ and t_f. The initial condition for the gradient flow is chosen as a piecewise constant function

where t_j is varied on a uniform grid between t₀ and t_f. A few examples of such initial conditions u(0, t) are shown as blue curves in Fig. 4(a), where t_j takes values uniformly in [3, 4], and all of the initial solutions converge to the same stationary solution represented by the red curve. In this example, we find six different stationary solutions {P_i: i = 1, …, 6}, each of which represents a different local minimizer of $Iσ$ given by (6); see Fig. 4(b). The paths {P₂, …, P₆} are the same up to translation by a multiple of the forcing period. P₁ is also a local minimizer of $Iσ$⁠, but its form depends on t₀, whereas the remaining optimal paths are insensitive to the choice of t₀ and t_f; see  Appendix C for details.

Figure 5 compares optimal paths, obtained from the gradient flow (11), for various values of σ. This includes setting σ = 0, which corresponds to using the FW functional. This example enforces boundary conditions so that all optimal transition paths start at the same initial point on $xl∗$ and terminate four periods later on $xu∗$⁠. We make the following observations:

1. Before crossing $xm∗$⁠, the optimal paths are concentrated along a similar trajectory, and they all leave $xl∗$ around the same phase of the forcing, marked by the green square. However, after crossing $xm∗$⁠, the optimal paths spread out in a way that is strongly dependent on the noise intensity.

2. The minimizer of I_FW stays near $xm∗$ for a full cycle before it exits to $xu∗$⁠. This is in contrast to the minimizers of $Iσ$⁠, where the path may transition quickly to the upper stable solution $xu∗$⁠.

3. Optimal paths for different σ cross $xm∗$ at different times and quickly relax to follow the deterministic flow after crossing.

4. For larger σ, a small gap between the optimal path and the stable periodic solutions before and after the transition is more visible. This points to a possible failure of the minimizers of the OM functional to adequately approximate the most probable path for larger σ since we expect the system to track the deterministic dynamics there. [We note that the contribution of I_OM2 to the functional (6) is more negative if the path is shifted from $xl∗$ to below $xl∗$ since f_x = 1 − 3x².]

The above observations can be understood by considering the relative scales of I_FW and I_OM2 both before and after crossing the unstable periodic solution for different regimes in Fig. 1, as described below.

Before crossing the unstable solution, the optimal path must go against the flow to reach that barrier, and, consequently, I_FW is necessarily large and bounded away from 0. For small σ, the contribution from the I_OM2 term is negligible in comparison with the I_FW term along this part of the optimal path. After crossing the unstable orbit, the optimal path can track the deterministic drift, and, hence, the contribution from I_FW can vanish for that portion of the path.¹⁹ As noted in Sec. II B, for the adiabatic regime, the I_OM2 term remains negligible throughout time and I_FW alone provides a good method for computing optimal paths.

In the nonadiabatic forcing regime, we consider, however, I_OM2 has the opportunity to play a role in shaping the portion of the optimal path after reaching the unstable orbit since contributions to I_FW may vanish there. This quality distinguishes optimal path dynamics in the adiabatic and nonadiabatic regimes. In particular, for the intermediate forcing regime, the I_OM2 term, which is related to the local expansion/contraction rate of the flow, selects when the optimal path departs the unstable orbit. Without this term, optimal paths computed using only I_FW (the blue curve in Fig. 5) track the unstable orbit for a full cycle under intermediate forcing. This is because the path can follow the deterministic dynamics without contributing to the value of I_FW. The inclusion of the I_OM2 term introduces a cost for noisy trajectories to remain near the unstable orbit, and this is important for computing optimal paths in our intermediate forcing regime. [Note that the contribution from I_OM2 is positive for portions of the path in the phase space band with $x∈(−1/3,1/3)$⁠.]

We numerically validate the path integral method by comparing optimal paths with the results of Monte Carlo simulations. The reported mean squared error associated with halving the number of samples quantifies our convergence test of the Monte Carlo results to the distribution. It is defined as

where $hN(m,k)$ is the fraction of samples in each bin m; see  Appendix D for more details about the convergence test and six comparisons for different values of (ɛ, σ).

Figure 6 illustrates results obtained with the path integral method compared to those from Monte Carlo simulations for parameters (ɛ, A, α, σ) = (0.1, 0.3, 0.15, 0.2). This figure shows the distribution of the sample paths on the (unwrapped) cylinder and two examples of distributions in x at t = −0.2 and 0, compared with the values of optimal paths (marked by red lines) at those times. These distributions are fairly concentrated as shown in the middle row of Fig. 6. We zoom in near the peaks of the distributions, in the bottom row, to better reveal the agreement between the optimal path and the Monte Carlo simulation results. Moreover, we notice that distributions broaden after t_m, the time when paths cross the unstable periodic solutions $xm∗$⁠, after which paths will simply follow the deterministic flow to $xu∗$⁠.

To quantify the discrepancy between the optimal path and Monte Carlo simulation results, we compute the mean absolute percentage error (MAPE) between the mode of time-evolved distributions and optimal paths for t ∈ [− 0.5, t_m], i.e., before crossing $xm∗$⁠. The MAPE is calculated as

where m_i denotes the mode of the distribution at each time step t_i, o_i is the value of the optimal path in x at t_i, $m¯=∑i=1kmi/k$ is the average of the mode of the distribution over t ∈ [− 0.5, t_m], and k represents the total number of time steps in t ∈ [− 0.5, t_m]. More comparison results for different values of ɛ and σ are summarized in Table I, and plots are presented in  Appendix D. Results from the Monte Carlo comparison suggest that the optimal path obtained from the OM functional gives a good estimate of the mode of the distribution prior to crossing $xm∗$ in this intermediate forcing regime. We also find that it does better than the FW functional alone; see comparison in Table I.

In this section, we use the optimal path method to investigate the intermediate forcing regime I of Fig. 1(a) and explore features of the deterministic problem that influence the tipping phase. Figure 5 demonstrated that transitions from the lower stable solution, $xl∗$⁠, to the unstable solution, $xm∗$⁠, may be insensitive to changes in σ, leading to a well-defined transition path associated with this part of the transition. Here, we explore the transitions between $xl∗$ and $xm∗$ more systematically for parameters in the nullcline regions R₂ and R₃, which were introduced in Fig. 3(b) for which the number of nullcline curves corresponds to 2 and 1, respectively.

Figure 7 summarizes results of a parameter study for α fixed at 0.1, σ ∈ [0.1, 0.4], and ɛ = 0.25 and 0.4. A is 0.4 in the left column and 0.7 in the right column, each associated with the different nullcline configurations that are in R₂ and R₃ of Fig. 3(b). We make the following observations:

1. For all four deterministic parameter sets, the optimal paths leave $xl∗$ near when a nullcline crosses this lower stable orbit, indicating a time when the flow changes direction from downwards to upwards near $xl∗$⁠.

2. The transition between $xl∗$ and $xm∗$ is concentrated around the same trajectory and varies little with the change of σ. Moreover, independent of noise intensity, the transitions from $xl∗$ to $xm∗$ all happen in the gap between the nullclines, which is where the deterministic flow is pointing upwards from $xl∗$⁠.

3. For ɛ = 0.25, the optimal paths arrive at $xm∗$ around the same time for different values of σ. However, for ɛ = 0.4, the times t_m differ by quite a bit. This is likely a result of the weaker deterministic flow associated with the larger value of ɛ.

4. Upon crossing the unstable periodic solution $xm∗$⁠, the optimal paths spread out and follow the deterministic flow to reach the upper stable periodic solution $xu∗$⁠.

We propose, based on the explorations summarized by Fig. 7, that a good deterministic predictor of the tipping phase t_l is the point where the flow near $xl∗$ changes direction from negative to positive. This point is associated with an appropriate nullcline crossing the lower stable orbit $xl∗$⁠. Moreover, for the flow geometries of regions R₂ and R₃, there is a gap between the nullclines; we conjecture that this gap serves as a passageway through which the most likely path squeezes. To further explore this claim, we focus on parameter sets in R₃ where there is a single nullcline, i.e., as in the cases of the right column in Fig. 7. We compute the fraction of the transition time from $xl∗$ to $xm∗$⁠, which lies above the nullcline (“outside the gap”), as a function of α, which controls the width of the gap; see Fig. 8. For α small (i.e., α near −0.135, the boundary between R₂ and R₃), the gap becomes very narrow, yet more than 80% of the transition from $xl∗$ to $xm∗$ occurs “inside the gap” below the nullcline, indicating the flow is upward. For α > −0.095 and below the saddle-node bifurcation, the entire transition from $xl∗$ to $xm∗$ is in this gap.

In this paper, we explored a notion of a “tipping phase” for a periodically forced Langevin equation (1) and extended the deterministic conditions that set it from the adiabatic to the nonadiabatic regime. Previous work in the adiabatic regime shows that the most likely transition occurs when the potential barrier height achieves its minimum,⁸ which for (2) is also the time when the separation between stable and unstable orbits is a minimum. Outside of the adiabatic regime, the two orbits exhibit roughly constant separation over time and the potential may lose its double well structure over an interval of time.

As a tool in our study, we used the path integral formulation to compute the optimal paths between two stable solutions. Specifically, we used the Onsager–Machlup functional whose minimizers we defined as optimal transition paths. The OM functional is the sum of the Freidlin–Wentzell rate functional and a second term that scales with σ²/ɛ, the latter of which controls the ratio of the noise to drift strength. We argued, and showed, that this second term in the OM functional can be significant in regime I in Fig. 1(a) when $ε=O(1)$ and σ ∈ (0, 1), and in particular, when the timescale ratio ɛ and noise intensity σ are comparable.

Our parameter study of Sec. III C suggests that the time when the flow direction changes sign, from negative to positive, near the starting lower stable solution is a robust estimator for the tipping phase, which we defined to be t_l, the time when the path leaves $xl∗$⁠. In some sense, this might be seen as an extension of the adiabatic regime results in which the transition happens at minimal barrier height. Specifically, in regions R₂ and R₃ where there are less than three nullclines, tipping starts when the double well structure first begins to break down, and optimal paths squeeze through regions where there is no barrier. More importantly, though, while in the adiabatic regime, tipping is nearly instantaneous (on the time scale of the forcing period) at the time point when the minimal potential barrier height is reached, in our case of an intermediate forcing regime, tipping takes much longer, making it necessary to define a tipping time more carefully. Our work suggests this estimator, t_l, and shows that long noisy paths attempt to make the transition by squeezing and maneuvering between portions of the nullclines.

We conclude with some future avenues of research that come out of this work and a discussion of some of the limitations of our approach that suggests open problems.

The Hamiltonian system (9) admits several periodic solutions (x^H(t), Ψ^H(t)) and has a rich bifurcation structure in σ and ɛ. For example, when σ = 0 and for typical parameters in our studies, there exist five periodic solutions of (9); three of these correspond to the solutions $x∗(t)$ of (2) and are given by $xH(t),ΨH(t)=x∗(t),0$⁠, while the other two are periodic solutions for which Ψ^H ≠ 0 (e.g., at A = 0, they satisfy $x=±1/3$⁠). We conducted a preliminary bifurcation analysis of (9) using AUTO-07P⁴⁰ and found that qualitative differences in the tipping trajectories in Figs. 7(d)–7(f) roughly align with saddle-node bifurcations of the Hamiltonian dynamics under changes of σ and ɛ. It would be interesting to explore whether the bifurcations in the Hamiltonian system can serve as a method for partitioning a (σ, ɛ) phase diagram into qualitatively different tipping behaviors. Moreover, near a deterministic bifurcation point, i.e., when periodic orbits $x∗$ collide, the deterministic dynamics enter a new scaling regime in which both noise and deterministic effects are comparable. In fact, near such a bifurcation point, the value of the integrand of I_OM2 may change sign. Consequently, we do not expect from our analysis that the intersection points of the nullclines with $xl∗$ will serve as a strong deterministic predictor for when tipping events will depart from $xl∗$⁠.

Our conclusions regarding the way the tipping phase aligns with a nullcline crossing are inherently tied to the low dimensionality of the problem we considered. Specifically, we considered a one degree of freedom, periodically forced system, for which an (extended) 1 + 1 dimensional phase plane analysis allowed us to gain insights. While our approach may not be limited to 1D systems, since we could base an investigation on the OM functional in higher dimensions, the natural analog of the nullcline is unclear. Is there a deterministic tipping phase indicator in some higher dimensional situations?

The benefit of the path integral approach is that it efficiently calculates optimal paths without performing Monte Carlo simulations. However, with this approach, we lose higher order information that is encoded in the distribution of tipping events. In contrast, in the Freidlin–Wentzell theory of large deviations, these quantities can be recovered from the so-called quasipotential in the limit of vanishing noise strength^21,23 or by asymptotics of the Fokker–Planck equation.¹⁹ A natural and open question is how to use the OM functional to obtain similar quantities in an asymptotic regime in which $ε=O(1)$ and noise is small (intermediate forcing rate) or ɛ ≪ 1 (adiabatic). This may be important in the context of the question we pose about a tipping phase for the following reason. If the distribution of tipping events is too broad, then there is no well-defined phase. It only gains meaning when the distribution is narrow compared to the period of the forcing. We expect that the distribution will narrow with decreases in σ, as illustrated in the contrasting histograms in Fig. 2, which were computed with different noise intensities.

Noise-induced transition has applications in a range of areas, such as eradication of infectious disease,^1–4,11 large fluctuations in chemical reactions,⁴¹ and climate systems,⁸ to name but a few. Our investigation arose, originally, in the setting of tipping points in a bistable conceptual model of Arctic sea ice, which has strong seasonal forcing. We wondered whether there is a season of the year when this important component of the climate system is most vulnerable to tipping from its current state of perennial ice cover to an alternative state where the Arctic is, year-round, ice-free.

Our preliminary investigations were based on an Arctic sea ice model due to Eisenman and Wettlaufer.⁷ Their model, an energy-balance column model of the Arctic, captures the strong seasonal variation in incoming solar and outgoing long-wave radiation over the course of the year and the key feedbacks associated with the ice-albedo and the sea ice thermodynamics. It is a 1D nonautonomous dynamical system, which shows bistability between the current Arctic state and the one where the Arctic is ice-free. We illustrate the noise-induced tipping phenomenon in this setting, using a version of the model due to Eisenman,⁶ with additive white noise. A stochastic version of this model has also been investigated by Moon and Wettlauffer,⁴² who consider both additive and multiplicative noise.

The Eisenman–Wettlaufer (EW) model⁷ has a prominent hysteresis loop under changes in greenhouse gases, akin to Fig. 3(a). Here, we consider the parameters in this model corresponding to the current state of the Arctic, which has perennial ice-cover and is bistable, in this model, with a perennially ice-free state; see the deterministic orbits in Figs. 9(a) and 9(b). Figure 9(c) presents distributions obtained from Monte Carlo simulations of the model with additive white noise. In contrast with our findings of the histograms for (1), summarized by Fig. 2, which showed the most prominent peak in the distribution for t_l, we find that the distribution of t_m is most peaked for the EW model, with a peak in April, which is, interestingly, well before the summer melting season. In Fig. 9(c), we shift all the tipping samples to cross the seasonal state in the third year and plot their mean (magenta curve). We note that, on average, the tipping events take many years to complete the transition, even with relatively strong noise intensity. Moreover, the noisy trajectories often stay in a neighborhood of the unstable seasonally ice-free state for almost a year or more before leaving for the perennially ice-free state, making it challenging to identify “a tipping phase,” despite the peaked distributions. These simulations appear to be in a fast forcing parameter regime where the intrinsic relaxation dynamics are slow compared to the yearly exogenous forcing time scale. We also note that the model has a natural nonsmooth limit, leading to a discontinuity boundary at E = 0, due to the transition from ice-covered to ice-free dynamics.⁴³ To further investigate its tipping events, considering the model switches that take place, with the seasonal effects, along the unstable seasonally ice-free state may be important.

We thank Julie Leifeld for some early discussions. A.V. thanks Björn Sandstede for his guidance in AUTO. We thank Paul Ritchie and Jan Sieber for introducing us to the path integral approach. M.S., Y.C., and A.V. are partially funded by the AMS Mathematics Research Communities [National Science Foundation (NSF) Grant No. 1321794]. The work of M.S. and Y.C. is funded in part by NSF Grant No. DMS-1517416. The work of A.V. and J.A.G. was supported in part by NSF-RTG Grant No. DMS-1148284, and A.V. is currently supported by the Mathematical Biosciences Institute and the NSF under Grant No. DMS-1440386.

Here, we provide a formal derivation of the Onsager–Machlup functional following the presentation in Chaichian and Demichev.²⁵ Let t₀, t_f ∈ ℝ with t_f > t₀ ≥ 0 and x₀ ∈ ℝ. For t ∈ [t₀, t_f] and g(x, t) a smooth function, consider the two continuous stochastic processes y(t), x(t) ∈ ℝ that satisfy the following stochastic differential equations interpreted in the Ito sense:

where σ > 0, y(t₀) = 0, x(t₀) = x₀, and W(t) is a Wiener process.

For a uniform discretization {t₀, t₁, …, t_N} of [t₀, t_f] of width Δt and realizations x(t) and y(t) of Eqs. (A1) and (A2), we define vector valued random variables x = {x₀, x₁, …x_N}, y = {y₀, y₁, …, y_N} by x_i = x(t_i) and y_i = y(t_i). It follows from the Kolmogorov extension theorem⁴⁴ that the probability density ζy for y is given by

In particular, if we let Δy_i > 0 be a sequence of spatial widths and $y¯$ a continuous function on [t₀, t_f] satisfying $y¯(t0)=0$⁠, then the probability realizations of Eq. (A1) pass through a sequence of “gates” $I=[y¯(t1)−Δy1,y¯(t1)+Δx1]×⋯×[y¯(tN)−ΔyN,y¯(tN)+ΔxN]$ is

The Wiener measure μ_y is defined as the unique measure on continuous functions satisfying (A4) for all temporal and spatial discretizations, i.e., all cylindrical sets.⁴⁴ The space of continuous functions is infinite dimensional, and, thus, a Lebesgue measure cannot be defined on it. However, as is customary in the physics literature, it is useful to think of the Wiener measure as heuristically having the following density with respect to Lebesgue measure:

where Z is the appropriate normalization factor.²⁵ Note that this heuristic notation results from considering the formal limiting density of Eq. (A3) as Δt → 0.

We can formally use Eqs. (A1) and (A2) to make a change of variables to construct the density ζx for x. In the integral form, it follows that

Using the composite trapezoidal rule, it follows that

Consequently, to the lowest order, the Jacobian matrix for this transformation is a lower triangular matrix with diagonal entries $1−Δt2gx(xi,ti)$ and thus to the lowest order

Since $exp⁡(−Δt2gx(xi,ti))=1−Δt2gx(xi,ti)+O(Δt2)$⁠, it follows, again to the lowest order, that

Finally, we obtain the following formal expression for the change of measure as Δt → 0:

Again, the representation of dμ_x as a Lebesgue measure is heuristic. However, using Girsanov’s theorem to rigorously compute the Radon–Nikodym derivative $dμxdμy$⁠, this result can be made precise for autonomous and nonautonomous drift.^45,46

To show that the optimal path is not sensitive to the choice of t₀ and t_f, we vary t₀ (t_f) by taking 16 uniformly spaced values in [0, 1] ([4, 5]) (Fig. 10). Initial guesses are in the form of (12) with the same intermediate time t_j = 2. In blow up box A (B), curves with different colors indicate that the paths start (end) at different positions on $xl∗$ (⁠$xu∗$⁠). In the second period, all paths converge to the same curve from $xl∗$ to $xu∗$⁠. After arriving at $xu∗$⁠, the paths start to separate shortly before the fourth period since we impose different final positions. Thus, given a large enough interval of time between t_l and t_u, optimal paths that start and end at different positions transition from $xl∗$ to $xm∗$ at the same forcing phase.

Here, we describe the procedure of Monte Carlo simulations. We use the Euler–Maruyama method⁴⁸ on the interval [t₀, t_f] with initial data $x0=xl∗(t0)$ to generate a sequence of samples X_j(t) for (1), each of which runs for K time steps (Fig. 11). We collect N samples that meet the following criteria for transitioning from $xl∗$ to $xu∗$⁠. Let $τj=∈ft{Xj(t)>xu∗(t)}$ denote the first time a sample path crosses $xu∗(t)$⁠. We define the “tipping events” to be {X_i} that satisfy {X_i} = {X_j:τ_j < t_f}. We then construct the distribution of the tipping events using M bins. The convergence test of the Monte Carlo simulations used to obtain distributions of the tipping events shown in Sec. III involves the following steps:

1. At each time step, t^(k) (k = 1, …, K), we divide the samples {X_i(t^(k)): i = 1, …, N} into M equal width bins in the x direction. We then record the fraction $hN(m,k)$ of the samples that fall into the mth bin. We pick M in the histogram based on the Freedman Diaconis rule,⁴⁹ and K so that Euler–Maruyama method converges.

2. We double the number of samples to 2N and repeat Step 1 to get $h2N(m,k)$⁠.

3. We then compute the mean squared error between $h2N(m,k)$ and $hN(m,k)$ (m = 1, …, M, k = 1, …, K), given by $E=∑k=1K∑m=1MhN(m,k)−h2N(m,k)2/(MK)$⁠.

4. If E < 10⁻⁵, we say that the Monte Carlo simulation has converged. If it has not converged, we double the number of sample paths and repeat the procedure.