Geometrical methods for analyzing the optimal management of tipping point dynamics

Natural resources are not infinitely resilient and should not be modeled as being such. Finitely resilient resources feature tipping points and history dependence. This paper provides a didactical discussion of mathematical methods that are needed to understand the optimal management of such resources: viscosity solutions of Hamilton–Jacobi–Bellman equations, the costate equation and the associated canonical equations, exact root counting, and geometrical methods to analyze the geometry of the invariant manifolds of the canonical equations.

As a general rule, I shall use "state" and "action" when developing theoretical relationships, and "pollution stock" and "loading" when discussing specifically pollution problems.
The differential equation x f x a x ẏ = ( , ), (0) = (1) describes the time evolution of the state, with y denoting its initial value. The time argument t of the variables x and a is suppressed, and x˙indicates differentiation of x with respect to time.
I assume particularly f x a a g x ( , ) = − ( ), that is, the natural decay g (x) depends only on the state level x. Two decay functions are illustrated in Figure 1. The dashed line corresponds to linear decay g x bx ( ) = .
( 3) Under time-independent actions a, the system settles at a steady state level x = a/b. Such a system is infinitely resilient, as the steady state level varies continuously with the action.
The "lake" decay (4) models resilience loss (Carpenter & Cottingham, 1997;Mäler, Xepapadeas, & de Zeeuw, 2003), for here the steady state is a discontinuous function of the loading. This is illustrated in Figure 2: for an initially unpolluted resource and small loadings, the steady state stock x increases almost linearly with a, as in Figure 2a. At the local maximum of the decay function, the system loses its resilience. If the loading increases past this value, the steady state jumps F I G U R E 1 Linear and lake decay functions discontinuously to a much higher value, as in Figure 2b. That figure also shows that decreasing the loading to its former value does not restore the pollution stock to its former level: it remains high.

| Economic valuations
For an ecologist, economic activities stress natural resources; for an economist, resources provide services. A lake provides freshwater, fishing, recreation, but it can also serve as a waste dump. To decide about how a natural resource should be best used, these services have to be valued. The cost flow of pollution is modeled as the shortfall cx − 2 of the income from an entirely unpolluted lake.
Phosphorus is a major pollutant of lakes, because it stimulates the growing of algae. It is also an important artificial fertilizer used in agriculture. In the economic lake literature, its benefits are modeled as a log .
Both lake and agriculture generate benefits. We assume a decision maker that regulates the use of artificial fertilizer, with the objective to maximize the discounted future benefits is associated to the stock level y. The dynamic programming approach consists in characterizing this value as the solution of a partial differential equation, the Hamilton-Jacobi-Bellman equation, and expressing optimal actions in terms of its derivatives. Making sufficiently strong assumptions, this equation can be derived readily. Assume first that the supremum in (6) is taken for an admissible action schedule a*(t) with corresponding state evolution x*(t).
For if equality did not hold, there would be an action schedule a(t) that is different from a* (t) on the time interval [s 1 , ∞) increasing the value. But that belies a*(t) maximizing J. Dividing (7) by s 1 − s, rearranging terms, and taking the limit s 1 → s yields For an arbitrary admissible state-action pair, Equations (7)-(9) turn into inequalities, as the realized benefit flow ℓ(x(t), a(t)) may be suboptimal. Hence the choice a = a*(s) maximizes the left-hand side of (9), which can be written as Note that calendar time nowhere enters the problem: dynamics and benefits only depend on state and actions, not explicitly on time. Hence the value is independent of time, V(s, y) = V( y), and after substituting x for y Equation (10) takes the form This is the stationary Hamilton-Jacobi-Bellman or HJB equation. It expresses, net of discounting, that for an optimal choice of a the decrease of future benefits equals the realization x a dt ( , ) ℓ of present benefits. Equation (11) holds at points where the value function is differentiable, if an optimizing action schedule exists. However, there are situations, important ones, where the value function is not differentiable. If we want the Hamilton-Jacobi-Bellman equation to describe these situations, we have to extend its meaning.

| THEORY OF HAMILTON-JACOBI-BELLMAN EQUATIONS
The compact way of writing the Hamilton-Jacobi-Bellman equation is For future use, note a particular feature of H. The expression in square brackets is a family, parametrized by x and a, of linear functions p x a pf x a ( , ) + ( , ).

↦ ℓ
The pointwise maximum of such a family is convex, therefore H(x, p) is convex in its second argument.

| Three difficulties
Working with Hamilton-Jacobi-Bellman equations present three technical difficulties.
The first difficulty we already have alluded to: the concept of the solution has to be enlarged to accommodate nondifferentiable value functions. Second, Hamilton-Jacobi-Bellman equations are implicit: the derivative V′ is not expressed explicitly as a function of x and V, which makes the equations awkward to work with. Finally, these equations are nonlinear and it is generally impossible to obtain closed form solutions.
For each difficulty, there is an appropriate mathematical tool. A nondifferentiable value function may be a viscosity solution of the Hamilton-Jacobi-Bellman equation. Crandall and Lions (1983), who developed this solution concept, showed, given suitable hypotheses, the value function to be the unique viscosity solution of the Hamilton-Jacobi-Bellman equation.
Implicit differential equations have been analyzed since the first days of calculus. The core idea is to take the derivative of the unknown function, the costate, as fundamental object. This works particularly well for the Hamilton-Jacobi-Bellman equation, as differentiation with respect to the state transforms it into a first-order quasi-linear equation of the costate. This equation may still have singularities, but these can be dealt with by a second transformation, resulting in the state-costate equations.
The remainder of this section discusses viscosity theory and the canonical equations. The next section analyses the lake management problem using Poincaré's geometrical methods.
3.2 | Viscosity theory 3.2.1 | The hedgehog problem A nondifferentiable value function already occurs if a hedgehog wants to get off a busy street. Let x ∈ [0, 1] be the hedgehog's position. If x = 0 or x = 1, it stands on the pavement and is safe from possible accidents. While 0 < x < 1, it incurs the possibility of being hit by a car. It prefers not to be hit, but it does not care about where or when: the expected benefit flow is negative and constant in space and time. We normalize it to −1. The hedgehog stands initially in position y, and it can move at a speed a: x a x ẏ = , (0) = .
Movement carries an exertion cost −a 2 /4. Let τ ≥ 0 be the first time the hedgehog reaches the pavement: subject to (13). The discount rate is set to zero. The hedgehog's Hamilton-Jacobi-Bellman equation reads as and it is complemented by the boundary conditions as the hedgehog incurs no damage on the pavement. Note the odd way the equation is written, with a minus sign in front of the function H: this will be significant later on.
Solving for the action a = 2V′(x) and substituting this into the Hamilton-Jacobi-Bellman equation reduces it to It is intuitively clear that the reasonable solution of this equation is which shows that the worst initial position is in the middle of the street. What is needed is a solution concept for the Hamilton-Jacobi-Bellman equation that results in this and only this solution.
Recall that an explicit first-order differential equation defines a direction field in the (x, V)-plane, the phase plane, as at every point the value of the derivative of V is specified uniquely through (17). In contrast, our Hamilton-Jacobi-Bellman Equation (15) defines a bi-direction field, as at every point (x, V) there are the two possibilities V′ = 1 and V′ = −1. Equation (15) has no everywhere differentiable solutions. For such a solution has to satisfy V′(x) = 1 or V′(x) = −1 for all x, but that is incompatible with the boundary conditions.
We might consider a function V to be a solution if it satisfies (15) almost everywhere. But then there are too many solutions, as any function satisfying the boundary conditions and whose derivative is 1 on a measurable set of measure 1 2 and −1 on its complement is a solution: for instance, the function W(x) = −V(x), with V as in (16), would be a solution.
The way out of this quandary is to apply an idea that originated in computational fluid dynamics. There numerical methods have great difficulties to compute the motion of nonviscous fluids, like water, as approximation errors generate spurious oscillations. The solution was to add a little "artificial viscosity," which suppressed the oscillations.
In dynamic optimization problems, a similar artifice is to add a little stochasticity to the motion: that is, we consider a slightly drunk hedgehog. Without delving here in the modeling details, this results in an additional second-order term to the Hamilton-Jacobi-Bellman equation. For the drunk hedgehog the equation reads as The equation can be integrated: V x x ε ( ) = − + 1 + log 1 + e 1 + e . Taking the limit ε ↓ 0 indeed selects our expected solution (16), as V ε (x) tends to V(x) uniformly in x. Denote by v( . Taking the zero-fluctuation limit of the stochastic problem yields the expected solution of the original Hamilton-Jacobi-Bellman equation. As the derivative of the stochastic value function converges everywhere to the derivative of the expected solution, except at x = 1 2 , Equation (18) implies that the second-order term ϵV x ( ) ″ ϵ tends to 0 everywhere except at x = 1 2 , and ( ) εV 1 ε ″ 1 2 → . That is, the limit V 0 (x) solves the Hamilton-Jacobi-Bellman equation almost everywhere, and where it fails to be differentiable, at x = 1 2 , it at least satisfies where v − and v + are the limits of V′ 0 as x tends to 1 2 from respectively below or above. In other words, ( ) v 1 2 is the derivative of a function that touches the graph of V(x) from below in x = 1 2 and for that value of

| Definition of viscosity solution
The example motivates the following definition. In contrast to the example, it is formulated in terms of the original equation, without recourse to an approximating higher order equation. The formulation is for functions on the real line, but it generalizes directly to functions of several variables.
if it is both a viscosity supersolution and a viscosity subsolution.
Clearly, if V is differentiable at x 0 , we can find test functions touching the graph of V at x 0 both from below and from above, implying that F(x 0 , V(x 0 ), V′(x 0 )) = 0.

| Relevant mathematical results
This section states the central mathematical results from the viscosity theory that will be used below in the analysis of the lake problem.
The first result states that value functions are viscosity solutions. The second states that supersolutions are larger than subsolutions everywhere. This then implies that viscosity solutions are unique and that the value function is the unique viscosity solution of the Hamilton-Jacobi-Bellman equation.
All these results hold under a number of technical assumptions. We formulate these for a situation where, unlike the hedgehog problem, there are no boundary conditions and the state space can be chosen a bounded open set. These assumptions are hypotheses for all the theorems in this section. In particular, they imply assumptions A 0 -A 4 in Chapter 3 of Bardi and Capuzzo-Dolcetta (2008), where also the proofs of the results can be found.
A1. The control set A is a compact subset of a finite dimensional vector space. A2. The space of admissible schedules A is the set of Lebesgue measurable functions a : are continuous, and hence bounded. A5. The functions ℓ and f satisfy a Lipschitz condition with respect to the first argument, that is, there is M > 0 such that for all x y , Ω ∈ and all a ∈ A The discount rate r satisfies r > 0.
As we are discussing only autonomous optimization problems, we consider the benefit functional and the evolution equation The Hamilton-Jacobi-Bellman equation associated to maximizing J subject to (20) is where H is given in (12).
is a viscosity solution of the Hamilton-Jacobi-Bellman Equation (21).
Theorem 2.2. Let u, v be bounded and uniformly continuous functions such that u is a viscosity subsolution and v a viscosity supersolution of (21). Then u ≤ v on Ω .
Compare Bardi and Capuzzo-Dolcetta (2008, Theorem 2.12). If there are two viscosity solutions u and v of the Hamilton-Jacobi-Bellman equation, then u ≤ v as u is a subsolution and v a supersolution, but also v ≤ u as v is a subsolution and u a supersolution. This implies u ≡ v, and the following theorem holds.
Theorem 2.3. The value function is the unique viscosity solution of the Hamilton-Jacobi-Bellman equation.
Hence, to find the value function of a dynamic optimization problem, it is sufficient to construct a viscosity solution of the associated Hamilton-Jacobi-Bellman equation. The construction will be effected by glueing together differentiable local solutions of the Hamilton-Jacobi-Bellman equation. In the hedgehog example we have seen that there are restrictions on how to glue local solutions together. The following theorem specifies these restrictions.
Theorem 2.4. Let x Ω ∈ , V(x) a viscosity solution of the Hamilton-Jacobi-Bellman equation, and φ: Compare Bardi and Capuzzo-Dolcetta (2008, Theorem 5.6). To apply the theorem, consider the situation p ↦ H(x, p) is strictly convex, and that the value function V(x), which is the unique viscosity solution of the Hamilton-Jacobi-Bellman equation, is continuously differentiable in a neighborhood of x, but that That is, the function V has a kink at x whose sharp point is pointing upward.
Then Theorem 4. implies, by continuity, that Hence if a value function is continuous and has piecewise continuous first derivatives, and the function H is strictly convex in its second argument, the value function satisfies the Hamilton-Jacobi-Bellman equations at points of differentiability and its derivative has to jump upwards at points of non-differentiability.
The following construction theorem states that these conditions are also sufficient.
Theorem 2.5. Let V : Ω  → be continuous and have piecewise continuous first derivatives. Assume that for every point x Ω ∈ , the following holds: Then V is a viscosity solution of the Hamilton-Jacobi-Bellman equation.
Proof. We have to show that V is a viscosity solution of the Hamilton-Jacobi-Bellman equation at points of nondifferentiability. Let x be such a point. There is a neighborhood of x such that V is differentiable for all points in the neighborhood, and

| Costate equation
We shall construct piecewise differentiable functions V that satisfy the Hamilton-Jacobi-Bellman equation at points of differentiability and the jump condition (22)  The solution strategy is to differentiate the Hamilton-Jacobi-Bellman equation with respect to the state variable, to obtain an equation that contains only the first and the second derivative of the value function, and that is linear in the second derivative, a so-called quasi-linear equation.
This equation can be analyzed with a dynamical system, the canonical system, that allows to construct the first derivative p = V′, or costate, of the value function. If the state space is higher dimensional, p is a vector-valued function whose components have to satisfy an integrability condition: this leads to the theory of symplectic geometry. In the present one-dimensional situation, integrability is automatically satisfied. To obtain the value function from the costate, the value has to be computed at one point in state space. Typically, this is done at steady states, where the objective functional can be evaluated directly.
Differentiating the Hamilton-Jacobi-Bellman equation yields This is an equation in p(x) = V′(x), the so-called costate equation By definition, the costate is the marginal increase in value if the state value is increased. For example, for a pollution problem, more pollution means less value, and the costate is typically a negative quantity.
As usual with first-order equations, it is convenient to rewrite this equation as a differential form and to note that any integral curve (x(s), p(s)) has to satisfy where α(x, p) is a real-valued function, only affecting the parametrization of the integral curves. For a dynamic optimization problem we have where a(p) is the optimal action for the given value of the costate that maximizes ℓ + pf on the right-hand side of (12). It follows that if α ≡ 1 on the right-hand side of (24), then dx ds H p x p f x a p dx dt = ( , ) = ( , ( )) = ∂ ∂ and consequently s = t + constant: the parameter s that parametrizes the integral curves is, up to an additive constant, equal to the time parameter. We find that trajectories of the canonical system are integral curves of the costate equation. Put differently, locally we can characterize the graph of the costate p(x) as an integral curve of the canonical system. This characterization will aid us in deriving properties of the costate function. Of course, the canonical system is the same as the state-costate system that features in the Pontryagin maximum principle. The above derivation shows how these arise naturally also in the "value function" approach.
The problem points of (23) are its singularities, which are the points (x, p) where p′ cannot be solved explicitly from the equation. These points are characterized by the equation Geometrically, these are all points where the tangent vector x p (˙,˙) to a trajectory is either zero or vertical, according to whether p˙is zero or not. We introduce this set formally as is positive below S x and negative above. If the tangent vector in a singular point is vertical, and if, when passing through S, the derivative x˙changes sign, then the trajectory "bends back" on itself and cannot represent the graph of a function. Such trajectories can, therefore, be removed from consideration. Therefore, the only possibility for the graph of the costate function p(x) to intersect the singular set S is in points where the tangent vector x p (˙,˙) is zero, that is, in equilibria of the canonical system. This section returns to the lake management problem and sets it up such that it is possible to apply the theory. To do so, some specific information is needed about the system response function that is summarized in the following proposition. Its proof is an exercise in calculus and is left to the reader.
1 is a local maximizer and x c 2 a local minimizer.
See Figures 2 and 4 for illustrations of decay functions with two critical points. The first and the third statement of the proposition bounds the parameter range of fragile systems, as it implies that there are no critical points if b > 3 3 /8. The second statement says that at the low critical point x c 1 , the system tips there toward the costly high pollution region, whereas at the high critical point x c 2 it tips back to the low pollution region. The last statement finally states that tipping back is impossible if the lake system is too fragile, that is, if b < 1 2 . Put briefly, there are no physical tipping points in the system if the system is robust, that is, if b > 3 3 /8, and there is no physical possibility from recovering from a bad tipping event if the system is very fragile, that is, if b < 1 2 . In the following, attention shall be restricted to the intermediate "reversible" case b > 1 2 : this is the situation where the lake can always be brought back to the unpolluted state x = 0 by reducing the inflow to 0.
The control set A is chosen such that both critical values of the system response are in the interior of A, and hence such that the lake can always be brought by management action either in the oligotrophic-nutrient-poor-region x x 0 < < c ℓ , or in the eutrophic-nutrient-richregion x x > h c . We therefore take A = [a ℓ , a h ] for constants 0 < a ℓ < a h such that a g x g x a < ( ) < ( ) < for all a ∈ A. Then no action can take the state to the complement of X (see Figure 4).
whose derivative φ(p) = Φ′(p) is given as Note for later use that the maximum in the definition of H is taken at a* = φ(p). In particular, if p ∈ C 1 ∪ C 2 , where C 1 and C 2 are the corner action regions, then the maximum is taken in a corner point of A, while if p ∈ I, the interior action region, it is taken in an interior point. The graphs of Φ and φ are shown in Figure 5. The Hamilton-Jacobi-Bellman equation reads as and the equation for the costate p(x) = V′(x) is φ p x g x p x r g x p x cx ( ( ( )) − ( )) ′( ) = ( + ′( )) ( ) + 2 .
The associated canonical system is When interpreting results, it is often convenient to write these equations in terms of the action variable a. This is of course only possible if p ∈ I, as a = φ(p) is invertible only there. For the lake, we obtain as state-action system valid for a(t) in the interior of the action space A.
In the following, we shall analyze this system for the fixed discount rate r = 0.03, which corresponds to a time horizon of log(2)/r ≈ 23 years, that is, approximately one generation. As this value has been widely used, the choice is made mainly to be consistent with the rest of the literature.

| Steady state analysis
By construction of the state space, all steady states of the state-action system are located in the interior of the product X × A of state space and action space. The steady states satisfy the two equations a g x a r g x cx = ( ) and or, equivalently, they are of the form (x, g(x)) with x a root of the equation Instead of solving this equation for x, which is hard, it is solved for c, which is easy, to obtain Plotting the graph of Ψ(x, b) for fixed values of b in a (c, x) diagram, as in Figure 6, shows for each value of c the number of steady states.

F I G U R E 5 Graphs of Φ(p) and φ(p)
The function Ψ not only determines the position of the steady states in the (c, x) diagram, it also determines their type.
Proposition 3.2. Let b > 1 2 and r > 0. If for a given parameter values (b, c, r) the statecostate system has a steady state x p (¯,¯) such that It can be shown that a stable steady state x a (¯,¯) of the state-action system corresponds to a weak local optimum of the benefit functional: that is, no loading schedule a(t) close to ā can achieve a better outcome. The proposition characterizes all such locally optimal loading schedules in terms of the function Ψ.
Proof. The derivative of Ψ reads as To determine the nature of the steady state, we use the state-action system (33), to obtain the derivative of the vector field at the steady state as The steady state Equations (34) can be used to replace r + g′(x) by 2cax and a by g(x), to yield Computing the determinant yields , then the determinant is negative, which implies that the eigenvalues of D are real and have opposite signs. This shows that the steady state is a saddle in this situation.
, the determinant is positive, and the eigenvalues have either both positive real part or both negative real part, so the steady state is either a repeller or an attractor. Now, the type is independent of the coordinates used to represent the system. We, therefore, pass to state-costate coordinates and consider the linearization of the statecostate system at the steady state. This is of the form As the trace of this matrix is r > 0-this is a property of general canonical systems-, we conclude that the sum of the eigenvalues is positive and the steady state is a repeller if > 0  Figure 6a). For c < 0.26, there is a single steady state x having values larger than 2. At the critical value, two additional steady states are generated in a saddle-node bifurcation at x = 0.60. For 0.26 < c < 1.09, three steady states exist. At c = 1.09 two steady states disappear in another saddle-node bifurcation at x = 1.09, and for larger values of c a single steady state remains around x ≈ 0.3.
If the resilience parameter b is increased, Figure 6b suggests that Ψ(x, b) is eventually monotonically decreasing in x. This means that the two saddle-node bifurcation points have disappeared in a cusp bifurcation, which necessarily occurs at a degenerate critical point.
To make these statements precise, we make use of the fact that Ψ(x, b) is a rational function, the quotient of two polynomials. In the following, we shall need information on the location of roots of polynomials that are derived from Ψ. There are exact algorithms to determine the number of roots of polynomials with rational coefficients in a given interval, which are far-reaching refinements of the well-known Descartes sign rule. They will allow us to formulate global results. These algorithms are too time-consuming to apply by hand, but they are implemented in symbolic computer algebra packages. Recall that Proposition 1. states that if and only if b < 3 3 /8 0.65 ≈ , there are multiple steady states under constant loading in the system, and hence physical tipping points. The present proposition is the economic analogue of this: it states that if and only if b < b cusp ≈ 0.77, there are multiple locally optimal steady states in the system. The apparent paradox that this can happen for values of b for which there are no physical tipping points is resolved by the remark that now time-varying pollution schedules are admitted.

Proof. A once degenerate critical value of x ↦ Ψ(x, b) is a cusp bifurcation value. That is, if a cusp bifurcation occurs for the parameter values
Note that since the denominator of Ψ(x, b) is nonzero for all (x, b) such that x > 0 and b > 1 2 , this property is shared with all derivatives of Ψ(x, b).
, with F and G polynomials in x and b, the latter two conditions in (38) , that is, the polynomial x ↦ F(x, b) has a double root x. Necessary and sufficient for this is that the discriminant of this polynomial has a root at b b =¯van der Waerden (1955, p. 93). For polynomials with rational coefficients, like in the present case, the number of roots in an arbitrary given interval can be determined exactly by a computer algebra package.
Appendix A shows the code of an implementation and its output. There are three results. First, the discriminant Δ(b) has a single root in the region b > 1/2. Second, the root is approximately located at b = 0.771419. Third, the interval [771/10 3 , 772/10 3 ] contains the root, which is the cusp bifurcation value b cusp . The first and third are exact results, which can in principle be replicated by hand.
The values for c cusp and x cusp can be determined similarly exactly, but that is not needed in the following. □ Of course, the interval containing the bifurcation value b cusp can be made as small as desired, and similar intervals can be found for the bifurcation values c cusp and x cusp .
From the theory of cusp bifurcations, it is known that two curves of saddle-node bifurcations in the (b, c) plane emerge from the cusp bifurcation point.
Proposition 3.4 (Saddle-node bifurcations). For b > b cusp , the system has a single steady state for all values of c > 0. For , such that the system has three steady states with positive values This proposition sharpens the result of Proposition 3. by giving tight restrictions on the set of parameter values for which there are multiple locally optimal states in the lake system. Geometrically, it has the shape of a wedge with vertex at (b cusp , c cusp ), and with the graph of c b ( ) Proof. Retaining the definitions of F and G from the proof of the previous proposition, by using a computer algebra package it is readily shown that F(x, 1) does not have a root in the interval x > 0. Hence the function Ψ(x, 1) is monotonically decreasing. If there were a value b b > cusp such that x b Ψ( ,ˆ) is not monotonically decreasing, by changing b from 1 to b, there is necessarily an intermediate has a double root for some x. This would be a degenerate critical point. But we know already that no such point exists for b > b cusp , hence b cannot exist and x ↦ Ψ(x, b) is monotonically decreasing whenever Similarly, in the region b b < <  , there are at most three interior steady states of the state-costate system (32), or, equivalently, of the state-action system (33), while for b ∈ (b cusp , ∞), there is only one.
This section constructs in the first situation solutions to the state-costate system that give rise to a viscosity solution of the Hamilton-Jacobi equation.
If the state-costate system (32) has a single saddle steady state, the Hamilton-Jacobi-Bellman equation has a differentiable solution, which is the stable manifold of the saddle. Recall that the stable manifold z x p = (¯,¯) is the set of all points z 0 such that the state-costate trajectory starting in z 0 tends to z as t → ∞. By the Invariant Manifold Theorem (Guckenheimer & Holmes, 1986;Hirsch, Pugh, & Shub, 1977), the stable manifold of a saddle in a twodimensional system is a smooth curve. This proposition is remarkable, as it gives a local criterium, in terms of the function Ψ, for the global solvability of dynamic optimization problem. The optimal management schedule steers the system from any initial state toward the steady state x.
Proof. By the hypothesis and Proposition 2., the unique steady state of the state-costate system x p (¯,¯) with p g x = −1/ (¯) is a saddle. From (37) we see that the vector (0, 1) cannot be an eigenvector of the matrix D. Then by the Invariant Manifold Theorem, locally around the steady state, the stable manifold of the saddle can be written as the graph of a function a w x =˜( ) s , or, equivalently, of a function p = w s (x). We extend the stable manifold by integrating the local manifold backwards in time.
The set Figure 7, is a backward trapping region: it is readily verified that on its boundary, all vectors are either tangent or outward pointing. Hence trajectories (x(t), a(t)) on the stable manifold can never leave B as t → −∞. In particular, the trajectories cannot cross the curve x a g ẋ = − ( ) = 0. Hence x(t) is monotonically increasing, if x t x ( ) <¯, or decreasing, if x t x ( ) >¯, and the trajectory can be represented as the graph of a curve a w x =˜( ) s for all x ∈ X. By construction, the function p(x) = w s (x) solves the costate Equation (31). The corresponding value function can be determined in two different manners. Either the function w s (x) is integrated, taking into account that at the steady state the integral J x a (¯,¯), where a w x = −1/ (¯) s , can be computed directly. Or it is determined directly from the Hamilton-Jacobi-Bellman Equation (30) as If the state-costate system has three steady states, then two of these have to be saddles. The stable manifolds of either of these saddles are then candidates for being the graph of the costate function. We, therefore, have to compare their values. There are three configurations of the relative positions of the two stable manifolds, one of which will be seen to be impossible on theoretical grounds. The first is that the stable manifold of one of the saddles is the graph of a function, that is, defined on the complete state space (Figure 8).
The second, which will be shown to be impossible, is that both stable manifolds define graphs of functions that are defined on the complete state space (Figure 9).
The last configuration is that locally both stable manifolds are described by functions defined on a subinterval of the state space; these subintervals may overlap ( Figure 10).
In the first configuration, where a stable manifold is equal to the graph of a function w X : i s  → , this gives rise to a smooth solution V(x) of the Hamilton-Jacobi-Bellman equation in the same manner as in the one steady state case: see Equation (39).
To see that the second configuration is impossible, we need the following auxiliary result. A region in the plane will denote an open set R whose boundary ∂ R is a piecewise differentiable curve. It is forward invariant under a dynamical system if every forward orbit with an initial point in the closure R of R is contained in R. is defined globally. Shown are the curves x˙= 0 (solid black), p˙= 0 (dashed), the stable manifolds W s 1 and W s 2 (gray), and the boundaries of the interior action region I (dotted) F I G U R E 9 Impossible configuration. Curves as in Figure 8 1 2 , r > 0 and c > 0 be fixed such that the state-costate equations have two saddle steady states x p (¯,¯) 1 1 and x p (¯,¯) 2 2 and one repeller x p (¯,¯) 3 3 . Let D 1 , D 2 ⊂ X be the closures relative to X of the maximal domains of definition of continuous functions w D : i s i  → such that the graph of w i s is a subset of the stable manifold W i s of the steady state x p (¯,¯) i i . Assume that neither D 1 nor D 2 is equal to X. Then the intersection D D D = 1 2 ∩ of D 1 and D 2 is nonempty.
Moreover, if D contains only a single point, then the Hamilton-Jacobi-Bellman equation has a smooth solution. If D is an interval of nonzero length, then the Hamilton-Jacobi-Bellman equation has a continuous but nondifferentiable viscosity solution, and D contains an indifference point.
This proposition characterizes the situation that there are two strongly optimal locally stable steady states. Here "strong" refers to the fact that all possible pollution schedules a(t) are admitted, not only those with values close to the steady state value. The indifference point separates two very different economic regimes: the one going to the low pollution steady state is characterized by large negative values of the costate p(x) = V′(x), hence a high incentive to reduce pollution through low loadings. The optimal management policy is characterized by a conservational approach and larger long-time benefits. The other economic regime is the one that steers the system to a high pollution steady state. That one is characterized by small negative values of the costate, and low incentives to contain pollution. The policy is an industryfirst approach that generates short-time benefits and larger long-time costs.
Proof. We have to investigate the various possible geometrical configurations of the stable manifolds. For the part of the stable manifold extending from the left-hand saddle point x p (¯,¯) 1 1 in the direction of increasing values of the state variable, these are illustrated in Figure 11.
Each orbit (x(t), p(t)) on this stable manifold locally parametrizes the manifold and accumulates on the saddle point for t → ∞. We shall classify the possible behavior of the orbit as t decreases.
By the Poincaré-Bendixson theorem, the distance of (x(t), p(t)) either increases to infinity as t → −∞, or it accumulates on a steady state, or on an invariant cycle.
x x s s x x 1 ′ 1 2↑ ↓ By Theorem 5. the function V X :  → defined as V(x) = V 1 (x) for x x x <≤ ℓ and V (x) = V 2 (x) for x x x < < h is a viscosity solution of the Hamilton-Jacobi-Bellman equation, and x is the indifference point. □

| Multiple actors
In this article, I have looked only at the single-actor situation. The main reason for this is that the multi-actor problem is not yet theoretically fully resolved. In Mäler et al. (2003) the situation that several decision makers affecting the pollution stock of a single lake by a pollution schedule conditioned on time-the so-called open loop information structure-has been analyzed. Largely, the analysis can be reduced to the single-actor situation, but with the cost, parameter c replaced by c/n, where n is the number of players. The more interesting feedback information structure, where the pollution schedule can be conditioned on the state, has been investigated in Dockner and Wagener (2014). There a "game" Hamilton-Jacobi-Bellman equation has been derived for the value function, which has again the same structure as the single-actor equation.
The implications for the dynamics are different, however. The jump in the costate value can occur even if there is only a single long-time steady state. Moreover, there are other, nonviscosity, solutions to the game Hamilton-Jacobi-Bellman equation that correspond to reasonable descriptions of the actors' behavior.