On the exact number of monotone solutions of a simplifed Budyko climate model and their different stability

We consider a simplified version of the Budyko diffusive energy balance climate model. We obtain the exact number of monotone stationary solutions of the associated discontinuous nonlinear elliptic with absorption. We show that the bifurcation curve, in terms of the solar constant parameter, is S-shaped. We prove the instability of the decreasing part and the stability of the increasing part of the bifurcation curve. In terms of the Budyko climate problem the above results lead to an important qualitative information which is far to be evident and which seems to be new in the mathematical literature on climate models. We prove that if the solar constant is represented by $\lambda \in (\lambda _{1},\lambda _{2}),$ for suitable $\lambda _{1}<\lambda _{2},$ then there are exactly two stationary solutions giving rise to a free boundary (i.e. generating two symmetric polar ice caps: North and South ones) and a third solution corresponding to a totally ice covered Earth. Moreover, we prove that the solution with smaller polar ice caps is stable and the one with bigger ice caps is unstable.


Introduction
The main goal of this paper is double: in a first step we study the exact number of solutions, depending on the parameter λ, of the discontinuous eigenvalue type problem where ω 2 is a given parameter and f (u) is the discontinuous function given by for some µ > 0, under the key assumption f 0 ∈ (0, 1), (1.2) with H(s) the Heaviside discontinuous function H(s) = 0 for s < 0, H(s) = 1 for s ≥ 0.
Here β is the maximal monotone graph of R 2 given by Notice that any solution of P (λ, f ) is also a weak solution of the multivalued problem P * (λ, β) −u (x) + ω 2 u ∈ λβ(u(x)) a.e. x ∈ (0, 1), The present paper can be considered as a natural continuation of the previous paper by the authors [3] in which the same program of research was devoted to the special case without any absorption term ω = 0. As we shall see, several of the sharp methods of proof used in [3] do not admit any easy adaptation to the apparently minor change introduced when assuming ω 2 > 0.
As mentioned in [3], Problem P (λ, f ) can be considered as a simplified version of some more general formulations arising in several different contexts: chemical reactors and porous media combustion, steady vortex rings in an ideal fluid, plasma studies, the primitive equations of the atmosphere in presence of vapor saturation, etc. We send the reader to the references collected in [3]. Nevertheless, our special motivation to study the case ω 2 > 0 was the consideration of problem P (λ, f ) as a simplified version of the so called Budyko diffusive energy balance models arising in climatology (see, e.g. [23], [18], [7], [25], [12] and a stochastic version in [11]). Although these models must be formulated on a Riemannian manifold without boundary representing the Earth atmosphere [12], the so called 1d-model corresponds to the case in which the surface temperature is assumed to depend only on the latitude component. The model considered in our previous paper [3] neglected the important feedback term arising when modeling the emitted terrestrial energy flux ( represented here, in a simplified form, by the term ω 2 u − f 0 ). In this way, we lead to a formulation similar to P (λ, f ) in which the spatial domain (0, 1) must be associated to a semisphere, the discontinuous function represents the co-albedo (with a discontinuity which is associated to the radical change of the co-albedo when the temperature is crossing −10 centigrade degrees: here represented by the value u = µ), the parameter λ the so-called solar constant, the boundary condition u (0) = 0 formulates the simplified assumption of symmetry between both semispheres and the condition u(1) = 0 represents the renormalized temperature at the North pole (i.e. we are assuming that u ≥ 0 in the rest of the hemisphere, and thus we must assume that µ > 0 although it represents −10 centigrade degrees).
We point out that the case in which the absorption term ω 2 u plays also an important role in the class of eigenvalue type problems P (λ, f ) associated to a discontinuous nonlinearity corresponds to the modelling considered by McKean [22] of the initial value problem for the FitzHugh-Nagumo equations which were introduced as a model for the conduction of electrical impulses in the nerve axon (see, e.g., Terman [26]).
We also recall that the very sharp bifurcation and stability results obtained trough the famous Crandall-Rabinowitz paper [4] requires in a fundamental way the differentiability of the nonlinear term. That was used in the very nice paper [18] to study the Sellers diffusive energy balance climate model in which β(u) is assumed to be at least a Lipschitz continuous function.
Results on the asymptotic behavior, when t → +∞, for the evolution energy balance model were obtained in [9] where it was also proved the general multiplicity of stationary solutions according the value of λ (see also [27] and [16] for other related results). A sharper bifurcation diagram, as a Sshaped curve was rigorously obtained in [1]. Nevertheless the method of proof in [1] uses the information obtained trough suitable zero-dimensional energy balance models and thus there is lacking of a more detailed information about the associated free boundaries generated by the solutions (given as the spatial points where T = −10). This kind of sharper information will be obtained here since no zero-dimensional energy balance model will be used in our proofs but only a direct analysis of the 1d-model.
One of the main difficulties to adapt the tools used in our previous paper [3] to the case in which ω 2 is not zero is the fact that the solutions of the ODE −u + ω 2 u = M may oscillate (in contrast with the case ω 2 = 0). As a matter of fact, there are several results in the literature indicating that the Budyko diffusive energy balance climate model admits an infinity of stationary solutions. That was shown in the [24] and [19] by considering a non-autonomous term λf (x, u) and in [13] for the mere autonomous case. The main goal of this paper concerns the study of non-oscillating solutions of the autonomous framework (which in the title is referred as "monotone solutions", as we shall explain below).
. Let u λ,µ be the monotone solution of P (λ, f ) given in iv).
Assume that u 0 − u λ,µ L ∞ is sufficiently small and let u(t, x : u 0 ) be the non-degenerate solution of P P * (λ, β, u 0 ). Then for some positive constant ε.
The proof of the instability of the stationary monotone solutions u λ,µ of the decreasing part of the bifurcation curve (λ, γ(λ)) can be obtained by different tools according the values of the parameter ω. A first possibility (as in [3]) is to study the sign of the first eigenvalue of the problem associated to the linearized equation where δ {x λ,µ } is the Dirac delta distribution at the free boundary point x λ,µ . When 0 < ω < 1 we can adapt the study made in [3] to prove that the principal eigenvalue η 1 of problem P η (x λ,µ : f 0 , λ) is negative.
Unfortunately, when ω ≥ 1 the above method of proof is not applicable and other arguments are needed. The following results use the sharp description of the equilibria given in Theorem 1.1 jointly to some monotonicity arguments.
In terms of the Budyko climate problem the above theorems lead to an important qualitative information which is far to be evident and which seems to be new in the mathematical literature on climate models. Let us extend, just by symmetry, the solutions of P (λ, f ) and P P * (λ, β, u 0 ) to the whole spatial interval (−1, 1) (this corresponds to consider the atmospheric surface temperature on the whole sphere instead on the south hemisphere). Assume the solar constant λ ∈ (λ 1 , λ 2 ). Then there are exactly two stationary solutions giving rise to a free boundary (i.e. generating two symmetric polar ice caps: North and South ones) and a third solution corresponding to a totally ice covered Earth. Moreover, among the two more realistic solutions (presenting two symmetric polar ice caps) the solution with smaller polar ice caps is stable and the one with bigger ice caps is unstable. Obviously, there are many aspects of the correct modeling of the energy balance climate models which were not taken into account in the simpler model analyzed in this paper. Nevertheless, it seems difficult to imagine that the stability properties of the more complex states which correspond to the monotone stationary solutions studied in our framework may be completely different to what is suggested here thanks to the "simplicity" of the model.
3 Stability of the increasing part of the bifurcation curve.
Proof of Proposition 3.1. For θ > 0, small enough, and let x λ+θ,µ the free boundary associated to the parameter λ + θ. By the continuity of the function (2.15) we know that there exists a small h(θ) > 0 such that x λ+θ,µ = x λ,µ + h(θ). Let us construct the upper barrier function ψ θ,δ in the following way: ψ θ,δ satisfies two different boundary value problems over different regions.

Instability of the decreasing part of the bifurcation curve
In this section, we shall prove Theorem 1.3 and Theorem 1.4 concerning the instability of lower branch u λ,µ .
Proof of Theorem 1.3 It is based on the fact that if the principal eigenvalue of P η (x µ,λ : f 0 , λ) is negative then the stationary monotone solution is unstable. The main idea of the proof is an adaptation of the similar result presented in [3] for the case ω = 0. If we consider the solutions v of the parabolic problem P P * (λ, β, v 0 ) with v 0 = u λ + φw 0 with w 0 smooth and φ small enough and to approximate v, as φ → 0, by functions of the form v(x, t) = u λ (x) + φw(x, t) with w(t, x) = e −νt U (x), with U solution of the eigenvalue problem P η (x λ,µ : f 0 , λ). Thus, using Proposition 4.1 in [3], we find that U satisfies , where ρ := η − ω 2 . We recall that in [3], we have studied the case ω = 0. In this case and when ρ = η := −τ 2 , we have showed that the free boundary x λ,µ ∈ (0, 1 2−f0 ) generate τ > 0 which gives a positive solution U.
Proof of Theorem 1.4. We shall use a very special construction of auxiliary initial data which leads, after a suitable convergence, to solution u λ,µ . Let θ > 0 small enough and take u 0 (x) := u λ−θ,µ (x).
Moreover, to check condition (1.4) observe that it is clear that u λ−θ,µ L ∞ (0,1) depends continuously with respect to θ. Indeed, this is exactly the continuity condition of the bifurcation curve given in Theorem 1.1. Notice, for instance, that the maximum point of u λ−θ,µ takes place at x = 0 and that, although at this point we only have a Neumann boundary condition, the fact that u λ−θ,µ verifies implies the above mentioned continuous dependence as a by-product of the continuous dependence of solutions of problem (4.21) with respect to the L ∞ (0, 1) norm of the right hand function and the continuous dependence with respect to the own interval of definition (recall that we know that the continuity of the function (2.15) implies that there exists a continuous function h(θ) > 0 such that x λ−θ,µ = x λ,µ − h(θ)).
An alternative proof of Lemma 4.1. Let v = ∂u ∂t , then v ∈ C 0 ((0, +∞) × (0, 1)) and verifies The maximum principle applies to (4.18) and we can conclude that v ≥ 0 (problems with measures of similar nature in [8] or [21]). Now, to conclude the proof of Theorem 1.4, we remark that (1.4) is already shown in Lemma 4.1.
Finally, if we define the sequence { u(t n ) L p (0,1) } tn>0 , ∀p > 1, then by Lemma 4.2, we see that this sequence is increasing and bounded. Hence, using a classical result we conclude that there exists ξ ∈ R, with ξ ≤ u λ,µ L p (0,1) such that u(t n ) u ∞ where u ∞ is a weak solution of the stationary problem corresponding to −u ∞ + ω 2 u ∞ = λf (u ∞ ) on (0, 1).
The following figure explains the dynamics of some initial data which are in small L ∞ (0, 1)−neighborhood of the instable solution u λ,µ : Figure 4: Dynamics of solutions corresponding to suitable initial data closed to the unstable equilibrium u λ,µ .

Remark 4.1
The proof of Theorem 1.4 also shows that the transient free boundaries corresponding to those initial data satisfy that x µ,λ (t) x λ,µ < 1 as t → +∞ and that x λ,µ (t) 0 as t → +∞.

Remark 4.2
The above proof generalizes and improves (with a different point of view) the results of [20] for the case of β a regular function.