Hopf bifurcation in a conceptual climate model with ice-albedo and precipitation-temperature feedbacks
Łukasz Płociniczak
math.DS physics.ao-ph

Spectral control for ecological stability
    Duccio Fanelli
    December 30, 2023
=========================================




	
In this paper we analyse a dynamical system based on the so-called KCG (Källén, Crafoord, Ghil) conceptual climate model. This model describes an evolution of the globally averaged temperature and the average extent of the ice sheets. In the nondimensional form the equations can be simplified facilitating the subsequent analysis. We consider the limiting case of a stationary snow line for which the phase plane can be completely analysed and the type of each critical point can be determined. One of them can exhibit the Hopf bifurcation for which existence we find sufficient conditions. Those, in turn, have a straightforward physical meaning and indicate that the model predicts internal oscillations of the climate. Using the typical real-world values of appearing parameters we conclude that the obtained results are in the same ballpark as the conditions on our planet during the quaternary ice ages. Our analysis is a rigorous justification of a generalization of some previous results by KCG and other authors. 


Keywords: climate dynamics, conceptual model, internal oscillations, Hopf bifurcation, dynamical system




§ INTRODUCTION

The story of our planet is a fascinating one having many spectacular moments and interesting phenomena. Focusing just on the climate, the data shows that the Earth has undergone a series of prolonged periods of low temperature during which the ice sheets descended to somewhat low latitudes - the ice ages. Apart some very extreme events - such as the hypothetical Snowball Earth (see for ex. <cit.>) - the last series of glaciations during the Quaternary period showed a quasi-periodic behaviour. These oscillations were characterized by a slow build-up of ice sheets and then a rapid deglaciation (see Fig. <ref>). This approximately regular recent behaviour of the climate is peculiar and they are several theories which try to explain it. The most popular one, founded by the research of Milutin Milankovitch <cit.> and later firmly established by the experimental findings of Hays, Imbrie, and Shackleton <cit.>, says that the astronomical variations in the Earth's orbit (axial tilt, eccentricity and precession) act as a pacemaker of the global temperature. 



Milankovitch theory is very attractive since the main frequencies of the astronomical forcing closely correspond to those observed in the geological data <cit.>. There are, however, certain problems with that simple explanation. First, the dominant period of the glacial oscillations during the last million years is about 100ka which closely corresponds to the period of orbital eccentricity changes. On the other hand, its can be easily calculated that the amplitude of that astronomical forcing is too small to drive the climate to the observed extent. Moreover, during the times before 1 Mya the dominant period of the temperature was approximately equal to 40ka. This peculiar bifurcation of the oscillation frequency cannot be ascribed to the Milankovitch theory (see <cit.>). As was suggested for example in <cit.>, the linear astronomical forcing can be treated as a tuning (or pacemaker) mechanism for a nonlinear internal oscillations of the climate. Some of the important ingredients of this mechanism are various positive and negative feedbacks from the ice sheets, precipitation, concentration of CO_2 and deep ocean circulation <cit.>. Moreover, the aforementioned change of oscillation period, the so-called Mid-Pleistocene Transition, can be explained by a particular bifurcation <cit.>. 

The various climate feedbacks can be analysed by the use of conceptual climate models which differ from the complex Global Climate Models (GCMs) by focusing only on some specific averaged properties of the Earth system <cit.>. These models are usually composed of several differential equations and can help to isolate the particular ingredients of the whole climate. To mention only a few examples, Saltzman and Maasch proposed and developed a physically based model describing glacial oscillations as a tectonically and astronomically forced limit cycle <cit.>. The main variables of the model are total ice mass, CO_2 concentration, and the North Atlantic deep water. Their theory has later been summarized in a book <cit.>. The mathematical treatment of this model has been recently conducted in <cit.>. Furthermore, a more elaborated box model of the ocean has been constructed in <cit.>. Recently, an interesting dynamical model exhibiting a subcritical Hopf bifurcation has been introduced in <cit.> modelling the ice volume, Antarctic ice sheet extent and the atmospheric concentration of CO_2. Moreover, in <cit.> the so-called minimal model for the ice-ages has been introduced as a forced van der Pol oscillator. Interestingly, some conceptual climate models exhibit mixed mode oscillations and this case has been rigorously investigated in <cit.>. 

Many approaches to modelling climate dynamics start from formulating the net energy balance averaged over the planet. These so-called Energy Balance Models (EBMs) have been initiated by seminal papers of Budyko <cit.> and Sellers <cit.> and nowadays are known by their name. Their various generalizations applied to many problems in climate dynamics have been investigated in a series of papers by North and his collaborators (for a summary see <cit.> and also <cit.>). There is also a recent and interesting rigorous study of the dynamics generated by Budyko-Sellers model supplemented by varying the ice line (extent of the ice sheet). Specific results can be found for example in <cit.>. Moreover, in<cit.> the ice energy balance notion has been supplemented by additional physical equations for the ice sheet movement and carbon cycle in an impressive dynamical system modelling these complex interactions. Finally, we mention the paper by Källén, Crafoord, Ghil <cit.> on which we base our reasoning. In that work the energy balance is coupled with ice sheet movement and a physically plausible limit cycle is numerically found. Further results on that model can be found in <cit.>.

In what follows we generalize the KCG model by introducing arbitrary functional forms of the ice-albedo and precipitation-temperature feedbacks. In reality, those functions cannot be determined explicitly by any measurement. We show that assuming only the physically obvious monotonicity and boundedness we are able to prove the results observed by previous Authors for the specific choices of aforementioned functions. We then consider a simplified version of the snow line model and rigorously justify that certain simplifications can be conducted in the nondimensional form of the system. Furthermore, we classify the type of each stationary point that can arise in this general situation. In particular, under conditions which we precisely state, one of them can exhibit a Hopf bifurcation yielding a limit cycle with an explicitly known first Lypaunov coefficient. In Section 2 we derive the model and rescale it to the nondimensional form. Section 3 contains the results of our analysis. 



§ MODEL DERIVATION



 §.§ Energy balance

We begin with a derivation of the governing equations which are statements of the energy and mass conservation. The initial step in the model formulation is a Budyko-Sellers type of energy conservation. In our work we assume the so-called zero-dimensional approximation meaning that we are concerned with globally averaged energy balance over sufficiently long time periods. Schematically, the energy equation has the form

    c dT/dt = q_i - q_o,

where T is the globally averaged temperature, c is the atmosphere heat capacity while q_i and q_o are, respectively, incoming shortwave and outgoing longwave radiative thermal fluxes. One can directly calculate that c defined in the above formula is approximately equal to m_a c_a / A_E, with m_a being the mass of the atmosphere, c_a specific heat of air and A_E Earth's area (see <cit.>).

Earth receives the energy in a high quality (mostly) shortwave solar radiation. Let Q denote the solar constant, i.e. the mean amount of Sun's irradiance per unit area of a plane perpendicular to solar rays. Assuming parallel ray approximation, Earth receives the amount of energy per unit time equal to π R_E^2 Q with R_E being the radius of Earth. Since 4π R_E^2 is the total Earth's surface area we have

    q_i = 1/4(1-α)Q,

where α is the mean terrestrial albedo (the fraction of reflected to absorbed radiation). It has to be noted that by no means should Q and α be regarded as spatially independent. For example, the time-average solar irradiance (so-called insolation) has a profound meridional distribution. This distribution can be calculated directly from the spherical geometry and has been done for example in <cit.>. One should point, however, that the horizontal distribution of heat by turbulent motion of the atmosphere takes place at a much faster time-scale than the one that interests us in our conceptual model (see <cit.>). Hence, it is not unreasonable to assume that the whole planet has a well-defined average uniform temperature. 

The aforementioned albedo α is a function of the underlying type of matter, for example ice reflects more light than a forest. In the present, when we globally average the albedo we obtain a value close to 0.3. Our foregoing considerations will be based on separating the continental α_c and oceanic α_o parts of the albedo

    α(T,l) = γα_c(l) + (1-γ) α_o(T),

where l is the latitudinal extent of the ice sheets and γ is a fraction of the area of continents to the whole surface of Earth. Of course 0≤α≤ 1. This model has been proposed in <cit.>. If more ice covers the land, the albedo should become larger. Therefore, α_c=α_c(l) will be later chosen to be an increasing function. On the other hand, the albedo of the ocean is usually taken to account for two extremes: the lack and the presence of sea ice. The formation of the sea ice is, in turn, strongly associated with the temperature: the colder the climate, the more sea ice is created increasing the oceanic albedo. For example, in <cit.> the albedo of the ocean has been taken to be a piecewise linear function, while in <cit.> a hyperbolic tangent smoothing approximation has been used. As we will show below, the exact form of the oceanic albedo is not relevant for the quantitative dynamics of the climate and hence we will use an arbitrary sigmoid function.

The outgoing longwave radiation (OLR) is Earth's main mean of energy loss. Because the atmosphere contains many absorbers sensitive to the long electromagnetic waves (hence the Greenhouse effect), the precise quantitative description of the OLR is a complex task (see for ex. <cit.>). Staying at the conceptual level we will adopt a semi-empirical approach in determining q_o. Since the magnitude of OLR closely agrees with incoming solar radiation, the planet is in a energetic equilibrium. We could thus use the Stefan-Boltzmann's law that states that the outgoing flux is equal to σ e^-Γ T^4. Here, σ is the Stefan-Boltzmann constant while Γ is an empirical parameter modelling the radiation absorption in the atmosphere - the so-called greenhouse parameter. This approach had been taken in a number of works, for example in <cit.>. We will, however, model the outgoing flux in a supposedly simpler way which dates back to Budyko's original work <cit.>. Mostly due to the short interval of relevant temperatures, say 250-320K, the sensible choice is to propose the following linear dependence 

    q_o = A + B T,

where A and B are constants determined from the data (to be given below in Tab. <ref>). It has also been previously argued that A can be related to the concentration of CO_2 in the atmosphere (see <cit.>). Finally, putting (<ref>) and (<ref>) in (<ref>) we obtain the equation of energy conservation 

    c dT/dt = Q/4(1-γα_c(l) - (1-γ)α_o(T))-A-B T,

where α_o is of the form (<ref>). Taking the l to be constant one can immediately see that the above equation enjoys multiple steady states which forms the basis of many fascinating phenomena associated with climate dynamics (see <cit.>). A careful examination of the corresponding bifurcation diagram shows hysteretic behaviour thanks to which the so-called tipping points of abrupt climate change are present <cit.>.



 §.§ Mass balance

By the mass balance we understand the growth and retreat of ice sheets under the influence of the climate variation. We will derive an equation describing the average behaviour of those great masses of ice. Our reasoning is based on the ice sheet model proposed in <cit.> (but see also <cit.>). This is a simplified model based on the plastic flow and isostasy. The more accurate and elaborated models based on Glen's law are surveyed in <cit.>. 



The overall picture of the ice sheet is presented on Fig. <ref>. The ice sheet is assumed to be zonally symmetric, i.e. have an essentially two-dimensional latitudinal profile. Moreover, we orient the x-axis pointing southward with origin at the ice sheet's half-width l. The northward rim x=-l is located at the Artic Ocean where the ice shelves form. Assume that the mass of ice of height h(x) above the sea level rests on the isostatically depressed bedrock. Following <cit.> we can assume that the rock density is three times the density of ice. Then, the buoyancy law states that the total ice thickness is equal to 3/2h(x). 

The ice sheet moves according to the plastic flow law, i.e. the whole mass behaves as a plastic fluid experiencing the yield stress τ_0. Since the normal stresses are hydrostatic, in the equilibrium we must have

    τ_0 = 3/2ρ_i g h |dh/dx|,

where ρ_i is the density of ice. This equation can be immediately integrated giving the parabolic profile

    h(x) = √(4τ_0 l/3ρ_i g(1-|x|/l))=H√(l)√(1-|x|/l),    H:=√(4τ_0/3ρ_i g).

The typical value of the yield stress and the corresponding height scale is given in Tab. <ref>.

The dynamics of the ice sheet is governed by an interplay of snow accumulation and melting. We assume that the northern part of the ice sheet is in equilibrium with its southern counterpart. Therefore, the mass balance of the whole sheet is governed by the ablation and accumulation in the region [0,l]. The amount of snowfall is closely associated with the temperature via the so-called snow line. The temperature falls with the height and thus there is a well-defined 0^∘ C isotherm under which the snow, if fallen, will melt. In the cited works this isotherm has been taken to be linear, i.e.

    h_iso(x) = h_0 + s(x+l),
	
where h_0 is the height at the Arctic ocean while s is the slope parameter. The northern height of the isotherm depends generally on the temperature or astronomical forcing. This has been taken into account in the literature. For example in <cit.> Authors correlated the variation in h_0 with the Milankovitch oscillations while in <cit.> a linear dependence on the global temperature has been prescribed.

The snowline position x=l_0 defining the part of the sheet which is nourished by the snowfall is defined as a solution of h(x) = h_iso(x). The mass balance can then be written as

    d/dt∫_0^l 3/2h(x) dx = a l_0 - m (l-l_0),

where a and m are respectively accumulation and melting rates. Now, if we assume that the ice sheet can grow or retreat, that is l=l(t), we can obtain

    d/dt∫_0^l 3/2h(x) dx = 3/2Hd/dt∫_0^l √(l)√(1-x/l) dx = 3/2H √(l) dl/dt.

Furthermore, the value of l_0 can be found by equating (<ref>) and (<ref>)

    h_0 + s(l_0+l) = H √(l)√(1-l_0/l),

which is a quadratic equation with the meaningful solution

    l_0 = H^2/s^2[-(h_0 s/H^2+s^2/H^2l+1/2)+√(h_0 s/H^2+2 s^2/H^2l+1/4)],

which together with (<ref>) gives the second dynamic equation

    3/2H √(l)dl/dt = (a+m) H^2/s^2[-(h_0 s/H^2+s^2/H^2l+1/2)+√(h_0 s/H^2+2 s^2/H^2l+1/4)] - ml   for  l_0 ≥ 0.

This is a rather complicated nonlinear equation but we will shortly see that an appropriate scaling and approximation will yield its accurate simplification. Also notice that nowhere in the above discussion we have mentioned the response of the bedrock to the evolution of the ice sheet. The moving mass will produce a delayed feedback from the lithosphere and hence, might provide an essential ingredient of the dynamics. This requires adding an additional equation to our system and a simplified version of it have been introduced in <cit.>. In this paper, however, we will consider only the instantaneous adjustment of the bedrock leaving the more general problem for the future work.

Two additions are necessary for the well-posedness of the model. First, it may happen that the whole southern side of the ice sheet is in the ablation zone. This can happen for sufficiently large sheets or by raising the 0^∘C isotherm. In this case the whole mass of ice becomes stagnant and our derivation is not valid (equation (<ref>) does not have a positive solution). Authors of <cit.> and <cit.> propose to assume that we can take l_0 = 0 and then, the mass balance can be written as

    3/2H √(l) dl/dt = - m l   for   l_0 < 0.

Another situation of (<ref>) failing to have a solution is when the negative snowline elevation h_0 becomes smaller than -sL. Then, the snow will accumulate on the ground and if l=0 the ice sheet will nucleate. For a sufficiently high snowfall rates a (or sufficiently small -h_0/s) the reasonable approximation is to assume (see <cit.>) that l = -h_0/2s and 

    3/2H √(l) dl/dt = - ah_0/2s  for   l<-h_0/2s,

while for l > -h_0/s the main equation (<ref>) still holds.



 §.§ Nondimensionalization

We are finally in position to nondimensionalize our model, i.e. equations (<ref>) and (<ref>). The summary of all of the physical quantities used in this model is given in Tab. <ref>, while a reference to all symbols appearing in the text is collected in Tab. <ref>.





First, we introduce the following scaling

    T (t) = T^* θ(τ),    l(t) = L^* λ(τ),    t = t^* τ,

where starred letters denote the appropriate scales while those in Greek correspond to the new dimensionless variables. The choice for T^* is straightforward. Since, as we mentioned before, the Earth is almost in thermal equilibrium where incoming solar radiation is balanced with OLR we can take

    T^* = Q/4 B,

which, using our typical data from Tab. <ref> is T^* = 195.55K. Then, the energy equation becomes

    c/B t^*dθ/dτ = 1+β-γα_c(λ)-(1-γ)α_o(θ) - θ,   β:=-4A/Q > 0,

where we have retained the same letters to denote (already dimensionless) albedo functions. Next, it is convenient to introduce the dimensionless snow line position in the following way

    l_0(t) = L^* λ(τ) λ_0(τ).

In this way λ_0 measures the position of the boundary of the ice sheet's accumulation level in terms of the total ice extent λ. We would also like to choose the ice sheet length scale L^* according to the 0^∘C isotherm. A quick look at (<ref>) lets us choose

    L^* = H^2/s^2,   ϵ = s h_0/H^2.

In this way we have

    λ_0(λ) = 1/λ[-(ϵ+λ+1/2)+√(ϵ+2λ+1/4)],

which implies the nondimensional form of the mass balance equation

    dλ/dτ = √(λ)((1+ξ)λ_0-1)   for  λ_0 ≥ 0.

where we have denoted the ratio of accumulation to ablation by ξ and chosen the ice sheet time scale

    t^* = 3/2H^2/m s.
 
Using the typical parameters from Tab. <ref> we can estimate that t^* = 33.2× 10^3 years, which roughly corresponds to the usual glacial oscillation time scale. Moreover, L^*=2.7× 10^4 km and ϵ = 0.1. Notice that ϵ can be thought as a small parameter. Having done all the scalings we can now formulate the final version of the temperature equation

    1/μdθ/dτ =  1+β-γα_c(λ)-(1-γ)α_o(θ) - θ,

where

    μ := 3/2B H^2/m s c,

is the free parameter measuring the ratio of ocean and ice sheet specific heats. We will also take this as our bifurcation parameter (what has also been done in <cit.>). This is because it is relatively easy to show that its variation produces a limit cycle via the Hopf bifurcation. And these in turn are of fundamental interest in climate science.

As for the modifications which take into accout l_0 < 0 and h_0 < 0 the nondimensionalization is straightforward. First, equation (<ref>) translates into

    dλ/dτ = -√(λ)  for  λ_0 < 0,

while (<ref>) becomes

    dλ/dτ = -ξ/2√(λ)ϵ  for  λ < -ϵ/2.
 
If λ≥ -ϵ/2 the main equation (<ref>) is valid.

To close the model we have to prescribe the albedos α_c,o and the melting to snowfall ratio ξ. The continental albedo should increase with the ice sheet extent. The simplest model of that is, of course, a linear dependence

    α_c(λ) = α_1 + α_2 λ,

where α_i are empirical constants. We have to keep in mind that in order of the above formula to produce a bounded albedo, the ice sheet extent has to be kept appropriately limited. The value α_1 is that of a clean ground while α_2 is determined from the continent's complete ice cover. If the ice is maximally developed, the continental albedo should become almost equal to that of the pure ice. Hence, we have the estimate

    α_2 ≥α_cmax-α_1,

where we have assumed that the maximal ice sheet extent is attained for λ=1. Thus, when we use the typical values from Tab. <ref> we obtain that α_1 ≥ 4. 

The oceanic albedo is a decreasing function of the temperature and, as we mentioned before, many authors proposed different versions of it: from stepwise functions to its continuous approximations. As we will show, the dynamics does not depend on the particular form of the albedo. 


	A function σ:ℝ→ℝ is called sigmoid if it is bounded and differentiable with a non-negative derivative. As a normalization one can take lim_x→±∞σ(x) = ± 1.
  
Therefore, it is reasonable to propose the following form of the oceanic albedo

    α_o(θ) = 1/2(α_+ +α_- + (α_+-α_-) σ_α(θ-θ_α/Δα)),   α_- ≥α_+,

where 0≤α_±≤ 1 are limits of α_o(θ) for θ→±∞, θ_α is the translation and Δα being the steepness parameter. The limiting case when Δα→ 0^+ corresponds to the discontinuous jump from α_- to α_+ at θ=θ_α. This form has been chosen for example in <cit.>.

In <cit.> Authors argue that ξ=ξ(θ) is an increasing function determined mostly by the snowfall rate. When the temperatures are low the oceans are mostly covered with ice. The evaporation is small, and hence the precipitation rate is low. Therefore, we can take

    ξ(θ) = 1/2(ξ_+ +ξ_- + (ξ_+-ξ_-) σ_ξ(θ-θ_ξ/Δξ)),   ξ_- ≤ξ_+,

where σ_ξ is an arbitrary sigmoid function and 0≤ξ_±≤ 1 are respective limits of ξ(θ) when θ→±∞. This concludes the derivation of our general model.



§ ANALYSIS



 §.§ Simplification

So far we have stayed at the full generality and did not make any simplifications corresponding to the particular configuration of the ice sheet and the temperature. It is crucial to notice that the λ-coordinate for the most interesting stationary point of the model is generally small. First, notice that by taking the value of h_0 estimated in <cit.> to be maximally of order of 1.2km, the reasonable magnitudes of τ_0 = 0.3 × 10^5 Pa and s = 0.4× 10^-3 we have ϵ≈ 0.1 which can be thought of as a small quantity. First, we will estimate the magnitude of the stationary value of the ice sheet extent. 

Notice that only the case of λ_0 ≥ 0 yields an interesting dynamics and therefore, in the following we will only focus on the analysis of (<ref>). It is relatively easy to obtain the estimates of the magnitude of the λ-nullcline on which the possible stationary points can lie. 


	Assume that
	
    -2+ξ/2ξ≤ϵ < 1/4ξ/2+ξ,
 
	where ξ=ξ(θ). Then the λ-nullcline of (<ref>) decomposes into two branches λ_1,2=λ_1,2(θ) for which 0≤λ_1 < λ_2. Moreover, the following holds for ϵ >0 
	
    (1+1/ξ)ϵ^2   ≤λ_1(θ) ≤ 4(1+1/ξ)ϵ^2, 
     2(1-1/2+ξ)ϵ   ≤ξ(1+ξ)/(2+ξ)^2-λ_2(θ) ≤ 3(1-1/2+ξ)ϵ,

	while for ϵ≤ 0 we have
	
    (1+1/ξ)ϵ^2+2(1+1/ξ)(1+2/ξ)ϵ^3    ≤λ_1(θ) ≤(1+1/ξ)ϵ^2,
    
    		-(1-1/2+ξ)ϵ   ≤λ_2(θ)-ξ(1+ξ)/(2+ξ)^2≤ -2(1-1/2+ξ)ϵ.



At the beginning assume that ϵ > 0. Observe that the λ-nullcline of (<ref>) with a help of (<ref>) can be determined from the equation

    1/λ[-(ϵ+λ+1/2)+√(ϵ+2λ+1/4)] = ζ,   ζ := 1/1+ξ.

Inspecting the left-hand side of the above we can notice that there exist exactly two solutions of the above if and only if the maximum of λ_0(λ) is greater than ζ. This maximum can easily be calculated to yield

    λ_0(λ_max) = 1-4ϵ/1+4ϵ  where  λ_max = 1/2ϵ(1+4ϵ).

Hence, two branches exist if and only if λ_0(λ_max) > ζ, that is to say ϵ < 1/41-ζ/1+ζ which is our assumption. 

Furthermore, since the equation defining λ_1 and λ_2 is essentially a quadratic, its solutions can be readily obtained in a closed, yet cluttered, form

    λ_1,2 = 1/2(1+ζ)^2[1-ζ-2(1+ζ)ϵ∓√((1-ζ)^2-4(1-ζ^2)ϵ)],

which after some manipulation can be written as

    λ_1,2 = 1-ζ/2(1+ζ)^2[1-21+ζ/1-ζϵ∓√(1-41+ζ/1-ζϵ)].

The next step is to estimate the square root in order to extract the leading order behaviour when ϵ→ 0^+. First, from the Taylor expansion we have

    √(1-41+ζ/1-ζϵ)≤ 1-21+ζ/1-ζϵ-2(1+ζ/1-ζ)^2ϵ^2.

From this we immediately obtain

    λ_1 ≥1-ζ/2(1+ζ)^2[1-21+ζ/1-ζϵ-1+21+ζ/1-ζϵ+2(1+ζ/1-ζ)^2ϵ^2] = ϵ^2/1-ζ.

Similarly, 

    λ_2 ≤1-ζ/2(1+ζ)^2[1-21+ζ/1-ζϵ+1-21+ζ/1-ζϵ]=1-ζ/(1+ζ)^2-2/1+ζϵ.

To find the other estimates notice that because of the convexity of the square root we have (or by some elementary manipulations)

    √(1-41+ζ/1-ζϵ)≥ 1-41+ζ/1-ζϵ.

Therefore,

    λ_2 ≥1-ζ/2(1+ζ)^2[2-61+ζ/1-ζϵ] = 1-ζ/(1+ζ)^2-3/1+ζϵ.

Finally, by the analysis of the auxiliary function 1-x -√(1-2x)-2x^2 we can derive the last essential inequality

    √(1-41+ζ/1-ζϵ)≥ 1-21+ζ/1-ζϵ-8(1+ζ/1-ζ)^2ϵ^2,

which gives us the last estimate for ϵ≥ 0. 

The case of ϵ < 0 is analogous. The maximum of λ_0 is always equal to 1 for ϵ < 0 and hence, both branches exist then and we can use the estimates obtained from the Taylor series. The first branch λ_1 can be estimated with a use of

    2(1+ζ/1-ζ)^2ϵ^2+4(1+ζ/1-ζ)^3ϵ^3≤ 1-21+ζ/1-ζϵ-√(1-41+ζ/1-ζϵ)≤ 2(1+ζ/1-ζ)^2ϵ^2,

for -1/21+ζ1-ζ≤ϵ≤ 0. The second solution λ_2 can be bounded by

    1≤√(1-41+ζ/1-ζϵ)≤ 1-21+ζ/1-ζϵ,   ϵ≤ 0.
 
This along with (<ref>) ends the proof.


Immediately we can see that 

    λ_1 = O(ϵ^2),   λ_2 = ξ(1+ξ)/(2+ξ)^2+O(ϵ)   as  ϵ→ 0,

and hence the branch converges to zero very fast with ϵ→ 0 while the larger is near a positive value. Henceforth, it seems that λ_2 is the more essential branch of the mass balance equation on which a possible critical point can lie. Moreover, because 0≤ξ≤ 1 we have

    ξ(1+ξ)/(2+ξ)^2≤ξ_+(1+ξ_+)/(2+ξ_+)^2≤2/9,

since the left-hand side function is increasing. Therefore, any critical point, if exists, has λ always close or smaller than 2/9 (for values from Tab. <ref> we have ξ_+=0.5 and thus the upper bound is 0.12). Because of that we can also consider it to be a small quantity. 

Now we are in position to make the final simplification of our governing equations. We want to utilize the fact that both ϵ and λ are small quantities. In order to simplify the expression for λ_0 we write

    λ_0(λ) = 1/λ[λ-(ϵ+2λ+1/2)+√(ϵ+2λ+1/4)].

Next, since ϵ+2λ is small we can expand the above square root into Taylor series to obtain

    λ_0(λ) = 1/λ[λ-(ϵ+2λ)^2+2(ϵ+2λ)^3+...]≈ 1-4ϵ-4λ-ϵ^2/λ,

where we have retained only two first terms of the above expansion. We can use this as an approximation of (<ref>) but we will make yet another simplification by taking ϵ→ 0. This is motivated by the above observation concerning the λ-nullcline branches λ_1,2. Finally, we arrive at our climate model which we repeat here for clarity. 

    {[ dθdτ =  μ[1+β-γ(α_1+α_2 λ)-(1-γ)α_o(θ) - θ]=:F(θ,λ),;               dλdτ = √(λ)((1+ξ(θ))(1-4λ)-1)=:G(θ,λ), ].
      for  θ> 0,    0 < λ≤1/4,

where α_o(θ), ξ(θ) are respectively given by (<ref>) and (<ref>) while all parameters μ, β, γ are positive. We will see that despite the apparent oversimplification of vanishing ϵ, this dynamical system possesses an interesting oscillatory behaviour which is in accord with the observations in paleoclimatology. Therefore, even the simplest coupling of realistic energy and mass equations produces self-sustained internal climate oscillations.

A note about validity of the model (<ref>) is in order. First of all, as ϵ→ 0 the smaller branch of the λ-nullcline disappears. The model should thus be valid only in the vicinity of the critical point lying on the λ_2 branch. On the other hand, the linear form of the continental albedo (<ref>) requires at most moderate variations of the ice sheet extent. Taking into account above arguments, we conclude that our simplifications are justified provided we do not allow for large excursions of λ.



 §.§ The phase plane

Since λ in (<ref>) appears only linearly in both of the given equations, we can immediately write the exact form of the nullclines (denoted by f and g)

    θ-nullcline:   λ =1/α_2[ 1/γ(1+β-(1-γ)α_o(θ)-θ)-α_1]=:f(θ),
       λ-nullcline:   λ = 1/4ξ(θ)/1+ξ(θ)=:g(θ).

An exemplary plot of the phase plane is depicted on Fig. <ref>. We can quickly compute the respective derivatives 

    f'(θ) = -1/γα_2((1-γ)α_o'(θ)+1),    g'(θ) = 1/4ξ'(θ)/(1+ξ(θ))^2.

Observe that g' is always positive while f' can change its sign depending on whether the oceanic albedo α_o has a sufficiently steep gradient. For the real-world parameters f has two local extrema and we denote them by θ_m and θ_M. We can see that f' is positive for θ∈(θ_m,θ_M) and negative otherwise.



Since α_o(θ) and ξ(θ) are both constructed with the sigmoid functions we can see that there can exist from one to five critical points of our system. With a simple geometrical reasoning we can distinguish several cases of their existence. 


	Let n be a number of critical points of (<ref>) and denote g_± := lim_θ→±∞ g(θ). If Δα is small enough for α_o'(θ) = (1-γ)^-1 to have two zeros θ_m and θ_M, then
	
		
  * n=1 if g_+ ≤ f(θ_m) or g_- ≥ f(θ_M) or (g(θ_m) < f(θ_m)  and  g(θ_M) > f(θ_M)).
		
  * n≥ 3 if f(θ_m) < g_-, g_+ < f(θ_M).
		
  * n=5 if and only if f(θ_m) < g_-, g_+ < f(θ_M) and there exists a point θ_c∈(θ_m,θ_M) such that f(θ_c)=g(θ_c) and f'(θ_c)<g'(θ_c).
	
	In the opposite case, we have n=1. 

Since the above can be justified by considering each case and plotting the nullclines, we omit the proof. Because of a high number of various parameters the result is stated in terms of the values of nullclines at readily computable points. This makes the presentation clear. Note also that in the above observation we include only the most meaningful situations. Of course, one can arrange the parameters in a way to produce an even number of critical points but this would require having the nullclines tangent at some points. In reality, this situation is highly improbable however, we will return to this case in Proposition <ref>. 

The (local) stability of the critical points can be categorized according to the gradients of f and g. Notice also that the creation and destruction of various stationary points depends on the parameters present in ξ and α_o (for ex. Δ T). Varying those can produce critical points through the saddle-node and pitchfork bifurcations. The precise conditions for them to occur are intrinsic to the specific nature of σ functions used in construction of ξ and α_o. We will treat all parameters except μ to be fixed and consider all possible bifurcations. 

First, we give a simple result concerning the linearisation.


	Let (θ_c,λ_c) be one of the stationary points of (<ref>) and by J denote its Jacobi linearisation matrix. Then
	
    J = (
    	μα_2 γ f'_c     -μα_2γ
    
    ξ_c√(λ_c)g'_c     -ξ_c√(λ_c),
    	)

	where we have used the short-hand notation ξ_c:=ξ(θ_c), f'_c=f'(θ_c) and g'_c:=g'(θ_c). Moreover, the eigenvalues of J can be written as
	
    r_1,2 = 1/2(μα_2γ f'_c-ξ_c/√(λ_c)±√((μα_2γ f'_c-ξ_c/√(λ_c))^2-4μα_2γξ_c/√(λ_c)(g'_c-f'_c))).



This is just a straightforward calculation. We start from collecting the formulas for partial derivatives of the vector field (F,G) defined in (<ref>)

    ∂ F/∂θ   = -((1-γ)α'_o(θ)-1),  ∂ F/∂λ = -μα_2γ, 
    ∂ G/∂θ   = √(λ)(1-4λ)ξ'(θ),  ∂ G/∂λ = 1/2√(λ)(ξ(θ)-12λ(1+ξ(θ))).

Since on the λ-nullcline (<ref>), in particular for the stationary point, we have λ=g(θ) the above partial derivatives of G at (θ_c,λ_c) reduce to

    ∂ G/∂θ = λ_c/√(λ_c)ξ'_c/1+ξ_c=ξ_c/√(λ_c)g'_c,  ∂ G/∂λ = 1/2√(λ_c)(ξ_c-3ξ_c/1+ξ_c(1+ξ_c))=-ξ_c/√(λ_c).

Moreover, the formula for the eigenvalues follows from the fact that

    r_1,2 = 1/2(trJ±√((trJ)^2-4 J)).

This ends the proof.

Notice that knowing (<ref>) we can easily find the direction of the vector field (F,G) since it vanishes on the nullclines and is monotone in their neighbourhood. 

Right now we are in position to subsequently classify the stability of all types of the stationary points that can arise in our system. At every critical point the nullclines f and g can cross transversalily or be tangent to each other. The transversality is the crucial property and can happen in two ways: either f'>g' or otherwise. This classifies the particular cases of stability.


	
	Let (θ_c,λ_c) be any critical point of (<ref>) for which f'(θ_c)<0. Then it is locally
	
		
  * a stable node for μ∈(0,μ_1] ∪ [μ_2,∞),
		
  * a stable focus for μ∈(μ_1,μ_2),
	
	where
	
    μ_1,2 := ξ_c/α_2γ(f'_c)^2 √(λ_c)(2g'_c-f'_c±√(g'_c(g'_c-f'_c)))>0.



Since f'_c<0 we immediately have tr J = μα_2 γ f_c' - ξ_c/√(λ_c) < 0 and thus the first term in both the eigenvalues (<ref>) is negative. Moreover, because g is always increasing we have f'_c < g'_c. Therefore, the quantity under the square root (<ref>) can become negative for a suitable choice of μ > 0. 

Immediately we notice that if the eigenvalues are conjugate complex the critical point is a stable focus since, as we mentioned before, Re  r_1,2 < 0. Next, observe that if the eigenvalues are real they have to be negative. The one with the minus sign is obviously smaller than zero while the other can be transformed into

    r_1 = 2μα_2γξ_c/√(λ_c)(g'_c-f'_c)/μα_2γ f'_c-ξ_c/√(λ_c)-√((μα_2γ f'_c-ξ_c/√(λ_c))^2-4μα_2γξ_c/√(λ_c)(g'_c-f'_c)) < 0.

Therefore (θ_c,λ_c) is a stable node. In any case, the critical point is locally stable.

We are left with finding the precise conditions for the transition from focus to the node to occur. According to (<ref>) we have to check whether the squared term chenges its sign. It is a quadratic equation for μ and can be written as

    (α_2γ f'_c)^2 μ^2 - 2 α_2γξ_c/√(λ_c)(2g'-f')μ + ξ_c^2/λ_c.

From the Vieta's formulas we immediately know that the zeros of the above are both positive. They can be easily calculated to yield (<ref>).


When we return to Fig. <ref> we notice that for the two outermost critical points we have f'_c<0 and hence, these are locally stable. The physical interpretation of this climate state is that we have either low temperature with a very small ice sheet or high temperature with a very large ice sheet. This situation is set up by the precipitation-temperature feedback and probably have to be corrected by allowing the snow line to vary with the temperature, i.e. ϵ=ϵ(θ). Since our simple model is tailored to investigate the other case we conclude that the two outermost critical points could not give a sufficient representation of the climate state. 

Now, we move to the most important and physical case where f is increasing. 


	Let (θ_c,λ_c) be a critical point of (<ref>) for which f'_c > 0. Then, the following holds.
	
		
  * When g'_c < f'_c the critical point is a saddle for all all μ>0.
		
  * When g'_c > f'_c the Hopf bifurcation takes place as μ passes
		
    μ_0 := ξ_c/α_2 γ f'_c √(λ_c).

		In this bifurcation a limit cycle is created with an approximate angular frequency 
		
    ω_0 = ξ_c/√(λ_c)√(g'_c/f'_c-1).

		Its stability can be determined from the sign of the first Lyapunov coefficient l_1(μ_0) which, for clarity of presentation, is given in the Appendix (bifurcation is supercritical for l_1(μ_0)<0 and subcritical otherwise).

		Moreover, the critical point locally is a 
		
			
  * stable node for μ∈(0,μ_1],
			
  * stable focus for μ∈(μ_1,μ_0),
			
  * unstable focus for μ∈(μ_0,μ_2),
			
  * unstable node for μ∈[μ_2,∞),
		
	where μ_1,2 are defined in (<ref>).
	


We start from the case of g'_c - f'_c<0. Then, from (<ref>) we can see that the term under the square root is always positive which implies that the eigenvalues are real for every μ > 0. Assume that μ is such that tr J ≤ 0, then the eigenvalue with a minus sign is negative while for the second we have a formula (<ref>) which shows that it has to be positive. Hence, the eigenvalues have opposite signs and the critical point is a saddle. The case tr J > 0 is analogous. 

Now, assume that g'_c - f'_c > 0 and notice that when determining the sign of the term under the square root in (<ref>) we arrive at exactly the same situation as before in Theorem <ref>. Hence, the critical point is a node for μ∈(0,μ_1] ∪ [μ_2,∞) and a focus for μ∈(μ_1,μ_2). The difference from this previous case is that its stability can change when μ increases. 

When μ=μ_0 the trace of the Jacobian matrix vanishes and the eigenvalues are purely imaginary. From the monotonicity of the quadratic we thus have μ_0 ∈ (μ_1,μ_2) and 

    Re  r_1,2(μ_0) = 0,   Im r_1,2(μ_0) = ± iω_0,

where ω_0 is defined in (<ref>). Moreover,

    .d/dμ|_μ=μ_0Re  r_1,2 = 1/2α_2 γ f_c' > 0,

hence the critical points cross the imaginary axis with nonzero speed. Therefore, by the Hopf bifurcation theorem <cit.> we conclude that a limit cycle is born when μ→μ_0. The stability μ-intervals can be obtained essentially in the same way as before with the use of (<ref>).

The last part of the theorem to prove is the form of the Lyapunov coefficient l_1(μ_0). First, we will transform (<ref>) into a standard form which helps us to utilize an explicit formula for l_1(μ_0). First, since at μ=μ_0 the Jacobian matrix has a conjugate pair of purely imaginary eigenvalues, we have

    J q = i ω_0 q,

where the eigenvector has a form

    q = (
    	[ 1/g'_c(1+i√(g'_cf'_c-1));                        1 ]).

Next, define the two transition matrices by

    S=(
    	[ q q ]),
    	  
    	T=1/√(2)(
    	[  1  i;  1 -i ]).

If we now introduce a new coordinate set by

    (
    	[ ψ; κ ]) = T^-1 S^-1(
    	[ θ-θ_c; λ-λ_c ])

then (<ref>) transforms into

    d/dt(
    	[ ψ; κ ]) =  
    	(
    	[    0 -ω_0;  ω_0    0 ])
    	(
    	[ ψ; κ ]) + 
    	(
    	[ P; Q ]),

where 

    (
    	[ P(ψ,κ); Q(ψ,κ) ])    = 
    	(
    	[                                                     G(θλ)+ω_0 κ; (g'_c/f'_c-1)^-1/2[G(θ,λ)- ξ_c/α_2γ√(λ_c)g'_c/f'_cF(θ,λ)]-ω_0 ψ ]) 
       where  (θ,λ)=(θ_c+(g'_c)^-1(ψ-κ(g'_c/f'_c-1)^1/2),λ_c+ψ).

Then, the first Lyapunov coefficient can be calculated via the formula (<cit.>, (3.4.29))

    l_1(μ_0)    = 1/8ω_0(P_ψψψ+P_ψκκ+Q_ψψκ+Q_κκκ) 
       +1/8ω_0^2[P_ψκ(P_ψψ+P_κκ) -Q_ψκ(Q_ψψ+Q_κκ)-P_ψψQ_ψψ+P_κκQ_κκ],

where subscripts denote partial derivatives evaluated at (ψ,κ)=(0,0). A straightforward calculation leads to (<ref>).

An exemplary plot of the limit cycle is presented on Fig. <ref> where we have used the parameters from Tab. <ref>. During our simulations we have noticed that although it is easy to control the position of the critical point, the stability of the limit cycle is much more difficult to anticipate. This is of course due to the rather complex form of (<ref>) which involves several competitive terms. However, due to the sigmoid nature of ξ and α_o their higher derivatives are very singular near their inflection point and thus, the presence of third derivatives in (<ref>) can dominate the behaviour of l_1(μ_0). The precise analysis of the Lyapunov coefficient is one of the objectives of our future work.



Looking at Fig. <ref> we immediately can see that although the position of the critical point is reasonable, the amplitude of the oscillations is smaller than the one predicted in the original work (<cit.>). During ice ages, the temperature oscillates through several degrees. The fundamental reason for such behaviour is putting ϵ = 0 at the beginning of our analysis. This forbids the snow line to follow the temperature changes and hence, to increase the temperature amplitude. This can be seen from simple physical reasoning. If the temperature is high, there is a larger snowfall at high latitudes promoting the growth of the ice sheet. Then, the large ice extent increases albedo which, in turn, lowers the temperature. In the case of ϵ≠ 0 the snow line gets lower and shrinks the ablation zone. Eventually, when the precipitation is low enough the melting dominates and ice retreats. Temperature rises forcing the snow line up and the oscillation continues so forth. On the other hand, when ϵ = 0 there is no additional mechanism which shortens the ablation zone and hence, the amplitude of oscillations is smaller than in the other case. 

On the other hand, when μ passes through μ_0 the amplitude of the oscillations grows until it reaches the size depicted on Fig. <ref>. This is in accord with the regime of model's validity - the excursions from the equilibrium should not be large. Moreover, if the parameters were tweaked for the nullclines to intersect precisely at one point, the limit cycle would not be attracted by another critical point. The oscillations would persist in a form of asymmetric relaxations-oscillations. The analysis of this case has been done in <cit.>.  

We stress, however, the fact that the period of the oscillations is of the right order that is present in the nature, that is 10-100k years (in our case - 40k years). Probably, it is possible to get more realistic results by varying the parameters (especially steepness Δα, Δξ and centres T_c, ξ_c) but our goal is to show that even for the most simple model of ice-albedo and precipitation-temperature feedbacks there is a stable limit cycle. 

Lastly, we state one more result about the possible type of the critical point. 



	Let (θ_c,λ_c) be a critical point of (<ref>) for which f'_c=g'_c>0. Then, if μ>μ_0 it is unstable while for μ<μ_0 its stability is the same as the stability of the ψ=0 point of the following differential equation
	
    dψ/dτ = 2 α _2 γμξ_c^2/(α _2 γ√(λ_c )μ  f'_c-ξ_c )^3(f_c”-g”_c)ψ^2 +  higher order terms.

	Therefore, if f”_c≠ g”_c the critical point is unstable. 


First by (<ref>), if f'_c=g'_c the determinant of the Jacobian matrix vanishes one of the eigenvalues is equal to zero. Then one eigenvalue vanishes while the other is equal to μα_2γ f'_c-ξ_c/√(λ_c). Now, if μ>μ_0 the latter is positive and hence, the critical point is unstable. 

Assume now that μ<μ_0. Then, one of the eigenvalues is positive while the other is zero. At this point we cannot infer about its stability and we will use the Central Manifold Theorem (see <cit.>). It simply states that on the phase plane there exist a graph of a function to which all sufficiently close trajectories converge exponentially fast. On this manifold, however, the dynamics is slower and governed by a smaller number of differential equations (in our case - one).

To find the center manifold first we have to diagonalize the system (<ref>) in order to find the direction in which the vector field has a zero eigenvalue. A simple calculation gives the eigenvectors

    p = (
    	[      -ξ_c √(λ_c ); -4 λ ^3/2ξ_c 'ξ_c ]),   q = (
    	[        γμα _1 f_c'; 4 λ_c ^3/2ξ_c 'ξ_c ]),

and the transition matrix

    S = ([ p q ]).

If introduce the following variables

    (
    	[ ψ; κ ]) = S^-1(
    		[ θ-θ_c; λ-λ_c ]),

the system (<ref>) becomes

    d/dt(
    	[ ψ; κ ])    = -ξ_c /4 λ_c  ξ_c  ξ_c '-4 α _2 γλ ^3/2μ  f' ξ '(
    	[ μ(4 λ_c ^3/2ξ_c '/ξ_c F(θ,λ)-α _1 γ  f'_c G(θ,λ));    4 λ_c ^3/2μξ_c '/ξ_c F(θ,λ)-ξ_c /√(λ_c )G(θ,λ) ]) =:
    	(
    	[ F(ψ,κ); G(ψ,κ) ]),
    	
       where  (θ,λ) = (θ_c+α _1 γμ  f'_cκ-ξ_c /√(λ_c )ψ, λ_c + 4 λ_c ^3/2ξ_c '/ξ_c (κ -ψ )).

From the Central Manifold Theorem there exist a function K=K(ψ,μ) which satisfies the following partial differential equation

    G(ψ,K(ψ,μ),μ) = K_ψ(ψ,μ) F(ψ,K(ψ,μ),μ),    K(0,0) = K_ψ(0,0) = 0.

Then, the dynamics of (<ref>) is captured by the equation

    dψ/dτ = F(ψ,K(ψ)).

Since we are dealing with the local behaviour near critical point we can expand the center manifold into Taylor series (this is the only meaningful coefficient)

    K(ψ) = c_2 ψ^2 +  higher order terms,

and plug it to (<ref>) to arrive at 

    c_2 = α _2 γ√(λ_c )μξ_c ' (4 λ_c^2 ξ_c ”-ξ_c^2 f_c”)+ξ_c ^2 ξ_c ”(ξ_c -α _1 γ√(λ_c )μ  f'_c)-8 λ_c  ξ_c ^2 (ξ_c ')^2/2 √(λ_c )ξ_c ' (ξ_c -α _1 γ√(λ_c )μ  f'_c)^2.

If we now plug the above expansion of the center manifold into (<ref>) we will arrive at

    dψ/dτ = α _2 γμ(ξ_c ^2 f'_c-4 λ_c ^2 ξ_c ')/ξ_c  (ξ_c -α _2 γ√(λ_c )μ  f'_c)ψ +  higher order terms.

However, by using (<ref>) and (<ref>) we have

    ξ_c ^2 f_c'-4 λ ^2 ξ_c ' = α_2 γμξ_c^2 (f'_c-ξ_c'/(1+ξ_c)^2) = α_2 γμξ_c^2 (f'_c-g'_c) = 0,

by our assumption. Therefore the leading order term in the above differential equation vanishes and we have to take into account higher orders of ψ. Using the above calculation and computing the next term in the Taylor series we obtain

    ψ' = 2 α _2 γμ(ξ_c ^3 f_c”+32 λ_c ^3 (ξ_c ')^2-4 λ_c ^2 ξ_c  ξ_c ”)/ξ_c (α _2 γ√(λ_c )μ  f'_c-ξ_c )^3ψ^2 +  higher order terms.

This can be further simplified by again using (<ref>) and noticing that

    4 λ_c ^2/ξ_c^2ξ_c ” -32 λ_c ^3/ξ_c^3(ξ_c ')^2 = 1/4ξ”_c/(1+ξ_c)^2 - 1/2(ξ'_c)^2/(1+ξ_c)^2 = g”_c.

From that there follows the assertion (<ref>). 




§ CONCLUSION AND FUTURE WORK

We have shown that even a very simple conceptual model describing ice-albedo and precipitation-temperature feedbacks can exhibit internal self-sustained oscillations. The model is by no means constructed to represent the actual state of the climate but rather to describe a simple mechanism which takes into account two important features of Earth system. 

Our results rigorously state that under some natural conditions, the temperature and the ice sheet extent oscillate at a period of O(10^4  years). The amplitude is, however, smaller than in reality but this behaviour can be anticipated from the very construction. Moreover, the results are completely independent from the particular choice of the functional form of the oceanic albedo and accumulation to ablation ratio (which are hard to determine experimentally). This generality strengthens the fact that the climate may indeed oscillate without any external forcing. 

Our future work will consists of considering the case of nonzero ϵ which would be a function of time and temperature in order to take into account the temperature-melting feedback of the ice sheets. The ultimate goal would be to augment the model with equations representing lithosphere response and CO_2 oscillations.



§ APPENDIX

The calculation of the first Lypanunov coefficient is a very technical point of the above considerations. We have to use (<ref>) and in order to do that, calculate the necessary second and third derivatives. For the function F we have

    ∂^2 F/∂θ^2   = α_2γ f”_c,   ∂^2 F/∂λ∂θ = ∂^2 F/∂λ^2 = 0, 
    ∂^3 F/∂θ^3   = α_2γ f”'_c,   ∂^3 F/∂θ^2 ∂λ = ∂^3 F/∂θ∂λ^2 = ∂^3 F/∂λ^3 = 0,

while the derivatives of G can be written as

    ∂^2 G/∂θ^2   = 4 λ_c^3/2/ξ_cξ_c”,   ∂^2 G/∂λ∂θ = 2 √(λ_c) (1-2 ξ_c )/ξ_c ξ_c',   ∂^2 G/∂λ^2 = -ξ_c/λ_c^3/2, 
    ∂^3 G/∂θ^3   = 4 λ_c ^3/2/ξ_cξ_c”'   ∂^3 G/∂θ^2 ∂λ = 2 √(λ) (1-2 ξ_c)/ξ_cξ_c”,   ∂^3 G/∂θ∂λ^2 = -(4 ξ_c +1)/√(λ_c)ξ_c ξ_c',   ∂^3 G/∂λ^3 = 3 ξ_c /4 λ_c ^5/2.

Having those we can compute the derivatives of P and Q in (<ref>). In order to facilitate this tedious process we have used the Mathematica scientific environment. The result is the following

    l_1(μ_0)    = 4 f_c^(3)λ_c ^2 ξ_c ^2+(f'_c)^2 (3 ξ_c ^2 g'_c-8 λ_c ^2 (4 ξ_c +1) ξ_c ')+8 λ_c ^3 (1-2 ξ_c ) f' ξ_c ”/32 λ_c ^2 ξ_c ^2 (f')^2 g'_c √(g'_c/f'_c-1)
       + λ ^2 ξ_c ^2 f_c”(4 λ_c ^2 ξ_c ”-ξ_c ^2 f_c”)+ξ_c ^2 (f_c')^3 (ξ_c ^2 g'+4 λ_c ^2 (2 ξ_c -1) ξ_c ')/8 λ_c ^2 ξ_c ^4 (f_c')^2 g_c' (f_c'-g_c') √(g_c'/f_c'-1)
       + 2 λ_c ^2 (f_c')^2 ((2 ξ_c -1) ξ_c ' (ξ_c ^2 g_c'+4 λ_c ^2 (2 ξ_c -1) ξ_c ')-2 λ_c  ξ_c ^2 ξ_c ”)-2 λ_c ^3 (2 ξ_c -1) f_c' ξ_c ' (ξ_c ^2 f_c”+4 λ_c ^2 ξ_c ”)/8 λ_c ^2 ξ_c ^4 (f_c')^2 g_c' (f_c'-g_c') √(g_c'/f_c'-1).

The sign of the above can now be numerically evaluated to infer the stability of a limit cycle. 




10
	
	Ash15
	Peter Ashwin and Peter Ditlevsen.
	The middle pleistocene transition as a generic bifurcation on a slow
	manifold.
	Climate dynamics, 45(9-10):2683–2695, 2015.
	
	Ber88
	André Berger.
	Milankovitch theory and climate.
	Reviews of geophysics, 26(4):624–657, 1988.
	
	Bud69
	Mikhail I Budyko.
	The effect of solar radiation variations on the climate of the earth.
	Tellus, 21(5):611–619, 1969.
	
	Cla06
	Peter U Clark, David Archer, David Pollard, Joel D Blum, Jose A Rial, Victor
	Brovkin, Alan C Mix, Nicklas G Pisias, and Martin Roy.
	The middle pleistocene transition: characteristics, mechanisms, and
	implications for long-term changes in atmospheric pco 2.
	Quaternary Science Reviews, 25(23):3150–3184, 2006.
	
	Cla02
	Martin Claussen, L Mysak, A Weaver, Michel Crucifix, Thierry Fichefet, M-F
	Loutre, Shlomo Weber, Joseph Alcamo, Vladimir Alexeev, André Berger,
	et al.
	Earth system models of intermediate complexity: closing the gap in
	the spectrum of climate system models.
	Climate dynamics, 18(7):579–586, 2002.
	
	Cru11
	Michel Crucifix.
	How can a glacial inception be predicted?
	The Holocene, 21(5):831–842, 2011.
	
	Cru12
	Michel Crucifix.
	Oscillators and relaxation phenomena in pleistocene climate theory.
	Phil. Trans. R. Soc. A, 370(1962):1140–1165, 2012.
	
	De13
	Bernard De Saedeleer, Michel Crucifix, and Sebastian Wieczorek.
	Is the astronomical forcing a reliable and unique pacemaker for
	climate? a conceptual model study.
	Climate Dynamics, 40(1-2):273–294, 2013.
	
	Dit18
	Peter D Ditlevsen and Peter Ashwin.
	Complex climate response to astronomical forcing: The
	middle-pleistocene transition in glacial cycles and changes in frequency
	locking.
	Frontiers in Physics, 6:62, 2018.
	
	Eng17
	Hans Engler, Hans G Kaper, Tasso J Kaper, and Theodore Vo.
	Modeling the dynamics of glacial cycles.
	arXiv preprint arXiv:1705.07387, 2017.
	
	Fow97
	AC Fowler.
	Glaciers and ice sheets.
	In The Mathematics of Models for Climatology and Environment,
	pages 301–336. Springer, 1997.
	
	Fow15
	AC Fowler.
	A simple thousand-year prognosis for oceanic and atmospheric carbon
	change.
	Pure and Applied Geophysics, 172(1):49–56, 2015.
	
	Fow13
	AC Fowler, REM Rickaby, and EW Wolff.
	Exploration of a simple model for ice ages.
	GEM-International Journal on Geomathematics, 4(2):227–297,
	2013.
	
	Fow11
	Andrew Fowler.
	Mathematical geoscience, volume 36.
	Springer Science & Business Media, 2011.
	
	Ghi81
	Michael Ghil and Hervé Le Treut.
	A climate model with cryodynamics and geodynamics.
	Journal of Geophysical Research: Oceans, 86(C6):5262–5270,
	1981.
	
	Ghi83
	Michael Ghil and John Tavantzis.
	Global hopf bifurcation in a simple climate model.
	SIAM Journal on Applied Mathematics, 43(5):1019–1041, 1983.
	
	Gil00
	Hezi Gildor and Eli Tziperman.
	Sea ice as the glacial cycles climate switch: Role of seasonal and
	orbital forcing.
	Paleoceanography, 15(6):605–615, 2000.
	
	Guc13
	John Guckenheimer and Philip Holmes.
	Nonlinear oscillations, dynamical systems, and bifurcations of
		vector fields, volume 42.
	Springer Science & Business Media, 2013.
	
	Hay76
	James D Hays, John Imbrie, and Nicholas J Shackleton.
	Variations in the earth's orbit: pacemaker of the ice ages.
	Science, 194(4270):1121–1132, 1976.
	
	Hog08
	Andrew McC Hogg.
	Glacial cycles and carbon dioxide: A conceptual model.
	Geophysical research letters, 35(1), 2008.
	
	Hyd00
	William T Hyde, Thomas J Crowley, Steven K Baum, and W Richard Peltier.
	Neoproterozoic snowball earthsimulations with a coupled
	climate/ice-sheet model.
	Nature, 405(6785):425–429, 2000.
	
	Kal97
	E Källén, C Crafoord, and M Ghil.
	Free oscillations in a climate model with ice-sheet dynamics.
	Journal of the Atmospheric Sciences, 36(12):2292–2303, 1979.
	
	Kuz13
	Yuri A Kuznetsov.
	Elements of applied bifurcation theory, volume 112.
	Springer Science & Business Media, 2013.
	
	Le83
	Hervé Le Treut and Michael Ghil.
	Orbital forcing, climatic interactions, and glaciation cycles.
	Journal of Geophysical Research: Oceans, 88(C9):5167–5190,
	1983.
	
	Len11
	Timothy M Lenton.
	Early warning of climate tipping points.
	Nature Climate Change, 1(4):201–209, 2011.
	
	Lis05
	Lorraine E Lisiecki and Maureen E Raymo.
	A pliocene-pleistocene stack of 57 globally distributed benthic
	δ18o records.
	Paleoceanography, 20(1), 2005.
	
	Maa90
	Kirk A Maasch and Barry Saltzman.
	A low-order dynamical model of global climatic variability over the
	full pleistocene.
	Journal of Geophysical Research: Atmospheres,
	95(D2):1955–1963, 1990.
	
	McG12
	Richard McGehee and Clarence Lehman.
	A paleoclimate model of ice-albedo feedback forced by variations in
	earth's orbit.
	SIAM Journal on Applied Dynamical Systems, 11(2):684–707,
	2012.
	
	McG14
	Richard McGehee and Esther Widiasih.
	A quadratic approximation to budyko's ice-albedo feedback model with
	ice line dynamics.
	SIAM Journal on Applied Dynamical Systems, 13(1):518–536,
	2014.
	
	Mil98
	Milutin Milankovič.
	Canon of insolation and the ice-age problem.
	Zavod za udžbenike i nastavna sredstva, 1998.
	
	Nad17
	Alice Nadeau and Richard McGehee.
	A simple formula for a planets mean annual insolation by latitude.
	Icarus, 291:46–50, 2017.
	
	Nor75
	Gerald R North.
	Theory of energy-balance climate models.
	Journal of the Atmospheric Sciences, 32(11):2033–2043, 1975.
	
	Nor81
	Gerald R North, Robert F Cahalan, and James A Coakley.
	Energy balance climate models.
	Reviews of Geophysics, 19(1):91–121, 1981.
	
	Nor83
	GR North, JG Mengel, and DA Short.
	Simple energy balance model resolving the seasons and the continents:
	Application to the astronomical theory of the ice ages.
	Journal of Geophysical Research: Oceans, 88(C11):6576–6586,
	1983.
	
	Oer80
	J Oerlemans.
	Model experiments on the 100,000-yr glacial cycle.
	Nature, 287(5781):430–432, 1980.
	
	Oer84
	Johannes Oerlemans and Cornelis J Van Der Veen.
	Ice sheets and climate, volume 21.
	Springer, 1984.
	
	Pai98
	Didier Paillard.
	The timing of pleistocene glaciations from a simple multiple-state
	climate model.
	Nature, 391(6665):378–381, 1998.
	
	Pai01
	Didier Paillard.
	Glacial cycles: toward a new paradigm.
	Reviews of Geophysics, 39(3):325–346, 2001.
	
	Pai04
	Didier Paillard and Frédéric Parrenin.
	The antarctic ice sheet and the triggering of deglaciations.
	Earth and Planetary Science Letters, 227(3):263–271, 2004.
	
	Pie10
	Raymond T Pierrehumbert.
	Principles of planetary climate.
	Cambridge University Press, 2010.
	
	Pie11
	RT Pierrehumbert, DS Abbot, A Voigt, and D Koll.
	Climate of the neoproterozoic.
	Annual Review of Earth and Planetary Sciences, 39, 2011.
	
	Plo19
	Ł Płociniczak.
	Asymptotic analysis of internal relaxation-oscillations in a
	conceptual climate model.
	arXiv preprint arXiv:1902.03467, 2019.
	
	Qui17
	Courtney Quinn, Jan Sieber, Anna S von der Heydt, and Timothy M Lenton.
	The mid-pleistocene transition induced by delayed feedback and
	bistability.
	arXiv preprint arXiv:1712.07614, 2017.
	
	Rob15
	Andrew Roberts, Esther Widiasih, Martin Wechselberger, and Christopher KRT
	Jones.
	Mixed mode oscillations in a conceptual climate model.
	Physica D: Nonlinear Phenomena, 292:70–83, 2015.
	
	Sal02
	Barry Saltzman.
	Dynamical paleoclimatology: generalized theory of global climate
		change, volume 80.
	Academic Press, 2002.
	
	Sal84
	Barry Saltzman, Anthony R Hansen, and Kirk A Maasch.
	The late quaternary glaciations as the response of a three-component
	feedback system to earth-orbital forcing.
	Journal of the Atmospheric Sciences, 41(23):3380–3389, 1984.
	
	Sel69
	William D Sellers.
	A global climatic model based on the energy balance of the
	earth-atmosphere system.
	Journal of Applied Meteorology, 8(3):392–400, 1969.
	
	Tzi06
	Eli Tziperman, Maureen E Raymo, Peter Huybers, and Carl Wunsch.
	Consequences of pacing the pleistocene 100 kyr ice ages by nonlinear
	phase locking to milankovitch forcing.
	Paleoceanography, 21(4), 2006.
	
	Van13
	Cornelis J Van der Veen.
	Fundamentals of glacier dynamics.
	CRC Press, 2013.
	
	Wal13
	James Walsh and Richard McGehee.
	Modeling climate dynamically.
	The College Mathematics Journal, 44(5):350–363, 2013.
	
	Wal16
	James Walsh, Esther Widiasih, Jonathan Hahn, and Richard McGehee.
	Periodic orbits for a discontinuous vector field arising from a
	conceptual model of glacial cycles.
	Nonlinearity, 29(6):1843, 2016.
	
	Wee76
	Johannes Weertman.
	Milankovitch solar radiation variations and ice age ice sheet sizes.
	Nature, 261:17–20, 1976.
	
	Wid13
	Esther R Widiasih.
	Dynamics of the budyko energy balance model.
	SIAM Journal on Applied Dynamical Systems, 12(4):2068–2092,
	2013.
	



