Normal forms for saddle-node bifurcations: Takens’ coefficient and applications in climate models

We show that a one-dimensional differential equation depending on a parameter μ with a saddle-node bifurcation at μ=0 can be modelled by an extended normal form y˙=ν(μ)−y2+a(μ)y3,where the functions ν and a are solutions to equations that can be written down explicitly. The equivalence to the original equations is a local differentiable conjugacy on the basins of attraction and repulsion of stationary points in the parameter region for which these exist, and is a differentiable conjugacy on the whole local interval otherwise. (Recall that in standard approaches local equivalence is topological rather than differentiable.) The value a(0) is Takens’ coefficient from normal form theory. The results explain the sense in which normal forms extend away from the bifurcation point and provide a new and more detailed characterization of the saddle-node bifurcation. The one-dimensional system can be derived from higher dimensional equations using centre manifold theory. We illustrate this using two examples from climate science and show how the functions ν and a can be determined analytically in some settings and numerically in others.

where the functions ν and a are solutions to equations that can be written down explicitly. The equivalence to the original equations is a local differentiable conjugacy on the basins of attraction and repulsion of stationary points in the parameter region for which these exist, and is a differentiable conjugacy on the whole local interval otherwise. (Recall that in standard approaches local equivalence is topological rather than differentiable.) The value a(0) is Takens' coefficient from normal form theory. The results explain the sense in which normal forms extend away from the bifurcation point and provide a new and more detailed characterization of the saddlenode bifurcation. The one-dimensional system can be derived from higher dimensional equations using centre manifold theory. We illustrate this using two examples from climate science and show how the functions ν and a can be determined analytically in some settings and numerically in others.
2022 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
The saddle-node bifurcation is one of the two generic codimension-one bifurcations of stationary points of differential equations (the other being a Hopf bifurcation). It occurs at parameters where the Jacobian matrix of a stationary point has a zero eigenvalue and two genericity conditions hold. Near a saddle-node bifurcation the local dynamics can be reduced to a one-dimensional differential equation on the associated extended centre manifold [1,2].
As parameters are varied to pass through the bifurcation, a pair of stationary points is created or destroyed. Restricted to the centre manifold, one stationary point is stable, while the other is unstable. On the centre manifold the system is locally topologically equivalent (in ways that will be made more explicit in §3) to the truncated normal forṁ y = ν − y 2 . (1.1) This truncated normal form is often introduced via coordinate transformations which push the other terms in the Taylor expansion of the family of systems beyond quadratic. These coordinate transformations can be given explicitly ( [2], section 3.3), yet the local equivalence is topological rather than differentiable except in the non-hyperbolic case when the eigenvalue of the Jacobian matrix is zero. In this case, successive coordinate changes can be made to remove nonlinear terms of all orders except quadratic and cubic [3], see also §2. This begs an obvious question: Why is it possible to obtain stronger results (differentiable rather than topological) in the apparently more exceptional non-hyperbolic case?
In this article, we show that the differentiable result extends beyond the bifurcation value provided a parameter-dependent cubic term is added to equation (1.1): (

1.2)
For this modified equation, the equivalence to the original equation restricted to the basin of attraction of the stable stationary point and, separately, to the basin of repulsion of the unstable stationary point is differentiable provided the parameter-dependent coefficients ν and a are chosen appropriately.
Here and later, we assume a reduction to the centre manifold has been made, so stability statements refer to stability restricted to the centre manifold. To be more precise, suppose the restriction to the centre manifold isẋ = f (x, μ), (1.3) and the saddle-node bifurcation occurs at x = 0 when μ = 0. Then f (0, 0) = f x (0, 0) = 0, and for genericity, we require f μ (0, 0) = 0 and f xx (0, 0) = 0. Without loss of generality, by changing of sign of the parameter and variable where necessary, we assume As shown in §5, to obtain a differentiable conjugacy between (1.2) and (1.3), it is necessary to have evaluated at x = μ = 0. This is exactly the quantity obtained for the non-hyperbolic case μ = 0 by Takens [3], see §2, so we refer to equation (1.5) as Takens' coefficient for the saddle-node bifurcation. It provides a quantitative measure of how far the system is from the truncated normal form (1.1) and breaks the time-reversing symmetry (y, t) → (−y, −t) of (1.1). It also reflects the degree of asymmetry in the two branches of stationary points in bifurcation diagrams, and a large value of a(0) indicates proximity to a cusp bifurcation. The remainder of this article is organized as follows. In §2, we revisit the formal power series approach to differentiable equivalence at μ = 0 and Takens' theorem for non-hyperbolic stationary points. Then, in §3, we bring together the main theoretical results on local equivalence of [4][5][6] in the form that will be needed in later sections. Thus, the results of § §2 and 3 are not new. In § §4 and 5, we use these results to modify the approach of [7] for bifurcations of maps to the continuous-time setting of saddle-node bifurcations. In §4, theorem 4.1 describes the smooth local equivalence in the case that no stationary points exist locally, proving that there is a C k conjugacy between the flow of f and the appropriate normal form. In §5, theorem 5.1 completes our analysis by showing that the conjugacy is C k on basins of attraction and repulsion of fixed points and collecting all the previous results for the saddle-node bifurcation into one statement. The new technical part of this result is obtained by using the implicit function theorem to match eigenvalues of the two stationary points to those of the extended normal form (1.2) and then using the classic differentiable conjugacy results of §3. This matching process determines the functions ν and a from f and its derivatives.
To illustrate the results, we have chosen two applications in climate modelling. In §6, we consider a simple temperature model of Fraedrich [8] and show how Takens' coefficient can be expressed in terms of physical attributes of the system. Then, in §7, we consider a version of Stommel's box model for ocean circulation [9]. This model is two-dimensional but has a slowfast time-scale structure, so perturbation techniques can be used to describe the leading order behaviour by a one-dimensional system for which Takens' coefficient can be evaluated explicitly. In §8, we explain how the reduction to a centre manifold can be performed, and in §9, we use this to evaluate Takens' coefficient for Stommel's model numerically. Finally, §10 provides concluding remarks.

Coordinate transformations: power series
In standard normal form analysis, the first step is to simplify when parameters are zero and then use a versal deformation argument to introduce general parameters to obtain local behaviour [1,10]. Coordinate changes are not performed away from the bifurcation point; instead topological equivalence arguments are used, often without detailed justification. On the other hand, to analyze saddle-node bifurcations using the implicit function theorem, the full strength of normal form theory is unnecessary and a shift of origin and scaling are used to remove the linear term in x from the Taylor expansion of the differential equation and then it is pointed out that the resulting system is locally topologically equivalent to the truncated normal form (1.1), see, e.g. [2].
In this section, we will recall the changes of coordinate at the bifurcation point for the saddlenode bifurcation leading to Takens' normal form theorem (theorem 2.1).

(a) Scaling the quadratic coefficient
With μ = 0, the general system (1.3) with (1.4) has the forṁ where here, and unless otherwise stated, derivatives are evaluated at x = μ = 0. To reach the modified normal form (1.2) with ν = 0, we first perform a linear change of coordinates, y = αx. We haveẏ (b) Persistence of the cubic term To the system (2.1), it is instructive to try and remove the cubic term via a subsequent coordinate change of the form: This inverts to y = z − βz 2 + O(z 3 ), sȯ and notice β is absent from the cubic term. Thus, the cubic term cannot be removed by this change of coordinates.
(c) Removal of quartic and higher order terms However, it is possible to remove higher order terms. Write (2.1) aṡ where k ≥ 4. Then with Since k > 3, we can choose which removes the z k term. This can be repeated for successively larger values of k removing all nonlinear terms of order four and above.

(d) Takens' theorem and general remarks
It is one thing to show that there is a formal power series which removes all terms higher than cubic, and it is quite another to show that this formal power series converges on a neighbourhood of the stationary point. Takens resolved this issue in the theorem quoted later. The aim of this article is to determine the correct formulation which allows us to accommodate μ = 0.
We are not aware of explicit C k versions of this theorem, but there are corresponding discretetime formulations [11], see also [7].

(e) Bifurcation theorems
In higher codimension problems, e.g. the Takens-Bogdanov bifurcation, and in some approaches to the Hopf bifurcation, a normal form is used at the bifurcation point and then a versal deformation argument is used to identify those small low-order terms (unfoldings) in parameter and phase space that imply all the topological behaviours close to the bifurcation point are realised [1,2,10]. For the simple codimension-one bifurcations, this approach is not necessary because topological equivalence is such a weak requirement. Presumably it is for this reason that theorem 2.1 is not often stated in the literature.

Equivalences and conjugacy
Formal proofs of even the local topological conjugacy results of bifurcation theorems are rarely given in textbooks, so in this section, we gather together the definitions and conjugacy results needed in the remainder of this article and state them in a form which simplifies their application in later sections. We emphasize that the results in this section are not new, and the novelty of this article lies in how these results determine properties of the saddle-node bifurcation described in § §4 and 5.

(a) Smooth conjugacies
This expresses the fact that y = h(x) is a change of coordinates witḣ Now let φ t (x) and ψ t (y) denote the flows induced byẋ = f (x) andẏ = g(y), respectively. An equivalent formulation of (3.1) is is a stationary point of g. By differentiating (3.1) and setting x = x * , we see that the two stationary points have the same stability coefficient, i.e.
This formulation is due to Sternberg ([6], theorem 4). In the general case of differential equations in R n , there are extra resonance conditions that need to hold between eigenvalues of the Jacobian matrices of f and g at the corresponding stationary points [4][5][6]. In the onedimensional case, there is no resonance, and more generally some differentiability is possible even with resonance if f is sufficiently smooth [12].
Two differential equations satisfying the C k -linearization theorem for the same value of λ = 0 are both conjugate toẏ = λy; thus, they are conjugate to each other. That is, we have the following result.

(c) Extension to basins of attraction and repulsion
We can now prove an adapted version of a theorem of Belitskii [13] for discrete-time dynamical systems that is key to the construction of the differentiable conjugacies for the bifurcations. A similar result is given in [7] for the discrete-time setting.
Let p ∈ (x j , u + ) and q = h(p) ∈ (y j , v + ). For each x ∈ (p, x j+1 ), there exists t x > 0 such that φ t x (p) = x and clearly t x is an increasing function of x with t x → ∞ as x → x j+1 . Define h(x) by using the two flows: from x evolve back to p, transfer to q = h(p), then evolve forward to h(x), i.e.
This is illustrated in figure 1. Since the flows are C k ( [14], theorem 0.8), by construction h(x) is C k and satisfies the conjugacy condition (3.2) for all x ∈ (p, x j+1 ). Thus, we have extended h to (u − , x j+1 ). The same argument in (x j−1 , x j ) allows us to extend h to (x j−1 , x j+1 ). (If f (x j ) < 0, then the argument is similar but with negative time.) Finally, we treat the case j = n (j = 1 is similar). As mentioned earlier, let p ∈ (x n , u + ) and q = h(p) ∈ (y n , v + ). As also mentioned earlier, we only treat the case f (x n ) > 0. The main difference here is that orbits ofẋ = f (x) reach the right endpoint of U in finite time, and similarly forẏ = g(y). To accommodate this, let t 1 , s 1 > 0 be such that φ t 1 (p) = u 1 and ψ s 1 (q) = v 1 . If t 1 ≤ s 1 , let x n+1 = u 1 and y n+1 = ψ t 1 (q) ≤ v 1 , and if t 1 > s 1 , let x n+1 = φ s 1 (p) ≤ u 1 and y n+1 = v 1 . Then the construction (3.3) produces a C k conjugacy for all x ∈ (p, x n+1 ), and as explained earlier, this is readily extended to (x n−1 , x n+1 ).
Note that the intervals (x j−1 , x j+1 ) are precisely the basins of attraction of x j if f (x j ) < 0 or the basin of repulsion of x j if f (x j ) > 0. Hence, another way to phrase theorem 3.3 is that there are local C k conjugacies on the basins of attraction and repulsion of corresponding stationary points.

Differentiable conjugacy and saddle-node bifurcations
The saddle-node bifurcation is associated with differential equations with a one-dimensional centre manifold reduction on which the equationẋ = f (x, μ) satisfies the conditions (1.4). In this section, we sketch the conditions that need to hold in order for C k conjugacies to exist, with proofs for the cases with no periodic orbits locally. Then, in §5, we apply these ideas to the normal form (1.2).
By the implicit function theorem, there exists a unique local branch of stationary points of the form This has two solutions in μ > 0 and none in μ < 0 for the choice of f μ > 0 and f xx < 0 made in (1.4). Now consider two C k families of differential equations, f parametrized by μ, and g parametrized by ν, both of which satisfy (1.4). Let φ t (x, μ) and ψ t (y, ν) denote their flows. If μ and ν are positive, then by theorem 3.3, there are local conjugacies if it is possible to choose ν = N(μ) such that the multipliers of the corresponding stationary points are equal. This is not a trivial condition. In §5, we use the implicit function theorem to prove that the coefficient a of the normal form (1.2) and the parameter ν can be chosen as functions of μ so that these conditions do hold, and so theorem 3.3 can be applied.
If instead μ = ν = 0, then both systems satisfy Takens' theorem (theorem 2.1) provided f and g are C ∞ . The following theorem accommodates the case in which μ and ν are negative. Theorem 4.1. Suppose that f and g are C k families of differential equations, f parametrized by μ, and g parametrzed by ν, both of which satisfy (1.4). There exist neighbourhoods U and V of zero and a C k function N(μ) with N(0) = 0 such that for all μ < 0 with |μ| sufficiently smallẋ = f (x, μ) on U is C k -conjugate tȯ y = g(y, N(μ)) on V.
Proof. If μ < 0 and ν < 0 with |μ| and |ν| sufficiently small, then the flows are decreasing in a neighbourhood of the origin. Given two such neighbourhoods U = (u 0 , u 1 ) and V = (v 0 , v 1 ), there exists T > 0 such that for all t > T there are functions μ(t) and ν(t) for which φ t (u 1 , μ(t)) = u 0 and ψ t (v 1 , ν(t)) = v 0 , for the first time. Moreover, since the flows are C k functions of the time and parameter (e.g. [15], sections 7.5 and 7.6 or [16], theorems 1.1.2 and 1.1.4), μ(t) and ν(t) are C k functions which tend to zero from below as t tends to infinity. Moreover, after further restricting the size of the neighbourhoods U and V and increasing T if necessary, standard comparison theorems (e.g. [15], section 2.7) imply that μ (t) and ν (t) are strictly positive. Hence, by the inverse function theorem [17], μ(t) can be inverted to obtain a C k correspondence ν = ν(t(μ)) = N(μ) between μ and ν.
At corresponding parameters, the systems are C k -conjugate on U and V by the C k flow box theorem ( [14], theorem 1.1). Equivalently, as this is trivial in one dimension, for all x ∈ U, there exists τ such that x = φ τ (u 1 , μ), and the C k conjugating function y = h(x) can be defined by N(μ)).

Normal forms and Takens' coefficient
In this section, we will prove the claims in §1. The analysis closely follows those in [7] for the discrete-time setting. There are four parts to the calculation: an asymptotic computation of the position of the stationary points of the general system (1.4), the equivalent computation for the normal form (1.2), the use of the implicit function theorem to identify functions ν(μ) and a(μ) for the parameter and coefficient of the normal form at which the corresponding stationary points have equal multipliers, and finally the use of theorem 3.3 to establish the existence of local differentiable conjugacies. Throughout this section, equation (1.4) can be used to write the differential equation for f (with all partial derivatives evaluated at x = μ = 0 and assuming that f is at least C 3 ) as follows: and by assumption, assuming f is at least C 4 . Differentiating (5.1) with respect to x and evaluating at x r gives the multipliers

(c) Equality of the multipliers
We require f (x r ) = g (y r ) for r = 1 and r = 2. These two equations are enough to determine ν and a as functions of μ. However, their leading order terms are equal (up to sign), and so to use the implicit function theorem, it is necessary to solve two related functions. Anticipating that n is proportional to m, define p by n = mp. Now (again a standard trick, see, e.g. [7,18]), let G r (a, p, m) = f (x r ) − g (y r ) and define Once again, from (5.4) and (5.7) with n = pm, and observe F 1 (a, p, 0) = 0 if p = p 0 and a = a 0 with With these choices of p 0 and a 0 F 1 (a 0 , p 0 , 0) = 0 and F 2 (a 0 , p 0 , 0) = 0.
To apply the implicit function theorem, the functions F 1 and F 2 must be at least C 1 ; hence, f must be at least C 4 . To obtain a unique local solution (valid in m > 0), the determinant of the matrix M of partial derivatives of F 1 and F 2 must not vanish. From (5.9) and (5.12),

(d) Local differentiable conjugacies
We can now put these calculations together with the strategy of §4 to determine the relation between the original map f and the normal form g (1. 2), which depends on the parameters ν and a.

y) on the basins of attraction and repulsion of the corresponding stationary points in U and V.
This simply collects together the results established earlier. Case (i) is theorem 4.1, note that if μ < 0, then any C k function could be used and a(μ) = a 0 would be a particularly simple choice. The condition a(0) = a 0 ensures continuity of the function with the solution to the implicit function theorem from case (iii). Case (ii) is Takens' theorem (theorem 2.1), and case (iii) follows from the results in this section together with theorem 3.3. The regularity of the parametrization follows from the fact that the implicit function theorem determining the functions ν and a is applied to a system of equations that is C k−3 . Note that this does not imply that these functions are C k−3 at the origin as the implicit function theorem is applied using the variable √ μ. We conjecture that theorem 5.1(ii) remains true in the C k case, k ≥ 4. This conjecture is certainly true in the discrete-time case [7,11].

(e) A remark on signs
Although f μ (0, 0) > 0 and f xx (0, 0) < 0 may be assumed without loss of generality for a generic saddle-node bifurcation, it is useful to have the constants p 0 and a 0 in a form which can be applied without making the explicit transformation to this case. Clearly regardless of the signs of f μ and f xx , but the sign of a 0 needs a little more thought. By considering the transformation x → −x, we see that in the differential equation, f xx → −f xx and f xxx → f xxx . Hence, a 0 is invariant under the transformation and no adjustment is necessary.

Fraedrich's temperature model
Saddle-node bifurcations are a mechanism for tipping points in climate dynamics. One such saddle-node bifurcation is exhibited by Fraedrich's model for temperature, T(t); see refs. [8,19] for details and justification. The temperature variation has a black body radiation term of order T 4 , a constant solar warming term, and a T 2 -term describing the variation of ice albedo with temperature. The one-dimensional model is 1 and the parameters are Here, the more fundamental constants are the Stefan-Boltzmann constant σ , the insolation I 0 , the oceanic thermal capacity c and the effective emissivity e SA . The constants a 2 > 1 and b 2 model the ice albedo effect: the albedo is a 2 − b 2 T 2 , and μ is a parameter that describes variations due to changes in the planetary orbit or insolation values. Reasonable values (taken from [19]) are The parameter μ is order one and acts as the bifurcation parameter.
so the solutions are Two solutions are created in a saddle-node bifurcation as μ increases through μ c = 4d/b 2 with temperature T c = 1 2 bμ c = 2d/b. Now let f (T, μ) denote the right-hand side of (6.1). A routine calculation shows that at the bifurcation point This means that the two quantities, p 0 and a 0 , that drive the normal form are related to the basic parameters by An interesting feature of this analysis is that Takens' coefficient, a 0 has two factors. The factor 1/a represents a scaling in time since in the new time τ = at the multiplicative factor in (6.1) is scaled to unity. The second (and indeed the last ratio in the expression for p 2 0 ) shows a nonlinear dependence on the level and sensitivity of the albedo effect to changes in global temperature.

Thermohaline circulation: perturbation theory
Stommel's two-box model of the oceanic circulation [9] remains a useful motivating example for the possibility of dramatic non-reversible changes in the climate due to global temperature rise. In the version used by Cessi [20] this can be written using a cubic nonlinearity instead of the modulus in the original model, so the non-dimensionalized temperature gradient x and salinity gradient y between the two boxes evolves according to the following equations: Here, α is a non-dimensionalized time constant, the ratio between the diffusive and temperature relaxation time scales, and is approximately 3600, m is the ratio of the diffusive time scale to the advective times scale and is approximately 7.5 and p is the non-dimensional freshwater flux with an average of order one. These estimates of magnitude are taken from [20], where p is allowed to oscillate stochastically about its mean and m is denoted by μ 2 , which we have changed to avoid confusion with general bifurcation parameters. Here, we take p to be constant and treat it as the bifurcation parameter (cf. [9]). Kuehn [21] described a slow-fast manifold approach to the same problem.
Since α is much larger than m and p, a perturbation theory approach can be taken with solutions expanded as power series in α −1 giving x = 1 + O(α −1 ) anḋ Writing y = y 0 + O(α −1 ), the perturbation equation for y 0 is obtained by ignoring the order α −1 terms, so, dropping the subscript 0 from here on, saddle-node bifurcations of (7.  The second equation of (7.3), which comes from setting the derivative of the right hand side of (7.2) equal to zero, is a quadratic for y = y(m). By solving it for y, we obtain and by substituting this into the first equation of (7.3), we obtain (after simplification) Note that p + < p − . If p ∈ (p + , p − ) and then there are three solutions, see figure 2. The two stable solutions correspond to those with highest and lowest salinity gradient. The ocean circulation currently lies on the solution with higher salinity gradient, and if p > p − , this is the only stable stationary point, lending weight to arguments that higher freshwater flux stabilizes the ocean circulation. However, if p decreases through p + , then this solution disappears and the solution would move rapidly to the stationary point with lower gradient and weaker circulation. This would have rapid and severe consequences for weather on the eastern European Atlantic coast. We therefore concentrate on the tipping point p + . If f represents the leading order term in the differential equation (7.2), then f p = 1, f yy = 4m − 6my, f yyy = −6m.
Observe f p > 0 and f yy < 0 at y + , and thus, the bifurcation occurs as p increases. The determining bifurcation numbers (evaluating at y + ) are

Centre manifolds
In higher dimensional problems, the one-dimensional system describing the saddle-node bifurcation is the projection of the flow onto a centre manifold. In this section, we use a centre manifold reduction to how the bifurcation numbers p 0 and a 0 can be computed for a two-dimensional system.