POSITIVITY OF SELF-SIMILAR SOLUTIONS OF DOUBLY NONLINEAR REACTION-DIFFUSION EQUATIONS

. The aim of this short note is to study self-similiar radially symmetric solutions of the scalar doubly nonlinear reaction-diﬀusion equation 1 on R n , where the parameters 1 < m,p,q < ∞ and 0 < λ < ∞ are ﬁxed. Particularly, for m < p < q < q c := p (1+ m − 1 n ) (where q c is Fujita’s critical exponent of blow-up) we show that there exist self-similar and radially symmetric solutions u , which do not blow up in ﬁnite time, but instantly become sign-changing for t > 0 inside some subdomain.

1. Introduction. Self-similar solutions are an important tool in the theory of partial differential equations, because they may capture certain aspects of the behaviour of solutions very well. Like Barenblatt solutions [2] for the porous medium equation such solutions are often source solutions, i.e. solutions with a Dirac measure at the origin as initial value. In this short note we study the self-similar and radially symmetric solutions u of the prototypical doubly nonlinear reaction-diffusion equation on R n for parameters 1 < m, p, q < ∞, 0 < λ < ∞, where u m−1 := |u| m−2 u, u q−1 := |u| q−2 u and (∇u) p−1 := |∇u| p−2 ∇u denote signed powers (and | · | denotes the absolute value resp. the Euclidean norm).
There are (at least) two different types of self-similar solutions: Solutions of Barenblatttype which exist for every time t > 0, and solutions which blow up in a finite time T > 0. In the pure diffusive case λ = 0 both types of self-similar solutions have extensively been studied in literature, see [18,Chapter 16] for a detailed analysis in the porous medium case p = 2, and [17,Part III] for the doubly nonlinear case. In presence of the reaction term λu q−1 , blow-up of solutions of (1) in finite time has been studied in the semilinear case m = 2 = p for q > 2 by Fujita in his seminal paper [5]. The semilinear case with a general reaction term f (u, ∇u) instead of λu q−1 has been discussed by [15]. A comprehensive overview of the porous medium case p = 2, m < 2, can be found in [14]. In this case, for a non-trivial non-negative bounded initial value blow-up in finite time occurs at almost every point for 818 JOCHEN MERKER AND ALEŠ MATAS m < q < 2 and at a set of measure zero for 2 < q < 2(1 + m−1 n ), while global solutions exist for q > 2(1 + m−1 n ) and sufficiently small initial values. The case m = 2, p > 2, of degenerate p-diffusion has been studied in [6]. In the doubly nonlinear case Fujita's critical exponent of blow-up, which in our notation is given by has been calculated in [7,Theorem 4.2]. The singular case 1 < p < m, p < 2, has been studied by [9]. For the case of an unbounded domain Ω = R n with noncompact Neumann boundary ∂Ω, see [1]. The case of a bounded domain Ω ⊂ R n under Dirichlet or Neumann boundary conditions is much more simple, there Fujita's critical exponent is identical with max(m, p), see e.g. [16] for the case m = 2, p > 2. Fujita's critical exponent and blow-up phenomena are surveyed in [10] and [3].
In this note we are particularly interested in the case of slow diffusion with a superlinear and superhomogeneous but subcritical reaction term, i.e. m < p < q < q c . We show that there exist self-similar and radially symmetric solutions which do not blow up in finite time but are global solutions vanishing at |x| = ∞. Further, provided that these solutions tend sufficiently fast to zero as |x| → ∞, they have a positive multiple of the Dirac measure at the origin as renormalized initial value. This seems to contradict some of the mentioned results in literature, e.g. [7,Theorem 4.2] states (at least in the doubly degenerate case 1 < m < 2 < p < ∞) that for m < q < q c the solution to every non-trivial non-negative initial value blows up in finite time. Further, for the case m = 2, p > 2, of p-diffusion in dimension n = 1 [6,Theorem 1] states that the solution to every non-trivial non-negative initial value blows up in finite time for 2 < q < p(1 + 1 n ). In fact, there is no contradiction, because the mentioned results assume at least implicitly that the non-trivial non-negative initial value is a bounded function. A suprising observation now is that although the Dirac measure Ψ := M δ 0 is positive for M > 0 in the sense of distributions, i.e. Ψ, φ ≥ 0 for every φ ∈ C ∞ c (R n ) with φ ≥ 0, the solution u(t) becomes instantly negative for every t > 0 at some subdomain of R n . Theorem 1.1. Let n ≥ 1, 1 < m, p < ∞, 0 < λ < ∞ and q c := p(1 + m−1 n ). If m < p < q < q c , then there exist a self-similar and radially symmetric solution u : (0, ∞) × R n → R of (1) such that u(t, ·) is sign-changing for every t > 0 and satisfies lim |x|→∞ u(t, x) = 0.
Moreover, if u(t, ·) µ−1 ∈ L 1 (R n ) for some t > 0 (and then in fact all t > 0), where 1 < µ < m is determined by q = p(1 + µ−1 n ), then u(t, ·) µ−1 converges to a multiple of the Dirac measure δ 0 at the origin in the sense of distributions as t 0.
The result that a solution u(t, x) to a positive measure as renormalized initial value becomes sign-changing for every t > 0, i.e. instantly negative at some points x ∈ R n , would have drastic consequences. Recall, that an evolution equation for u is said to preserve positivity, if u 0 ≥ 0 for the initial value at t = 0 implies u(t) ≥ 0 for every t > 0, and that the weak comparison principle is said to hold, if u 0 ≤ũ 0 implies u(t) ≤ũ(t) for every t > 0. Thus, preservation of positivity and the weak comparison principle would be violated, if a solution u(t, x) to a positive measure as renormalized initial value becomes sign-changing for every t > 0, i.e. for general measure as renormalized initial value equation (1) would neither preserve positivity nor admit a weak comparison principle.
In contrast, for initial values which are not so irregular as above (where u µ−1 0 = δ 0 for µ < m and particularly the power u m−1 0 does not have a rigorous meaning), but more regular in the sense that u m−1 0 ∈ L 1 (R n ), positivity is preserved and a weak comparison principle holds. In fact, under the same condition q > m as in Theorem 1.1 the right hand side f : v → λv q−1 m−1 of (1) is locally Lipschitz continuous as a function of v := u m−1 , thus for solutions u m−1 ,ũ m−1 ∈ C(0, T ; L 1 (R n )) to initial values u m−1 0 ,ũ m−1 0 ∈ L 1 (R n ) which a priori satisfy an L ∞ -bound |u(t, x)|, |ũ(t, x)| ≤ M for a.e. t ∈ (0, T ), x ∈ R n , the L 1 -contraction principle of [13] implies Outline. In section 2 we review Barenblatt solutions of doubly nonlinear diffusion equations, i.e. of (1) with λ = 0. In the main section 3 we calculate the global self-similar and radially symmetric solutions u of (1), derive their properties and prove Theorem 1.1.

Barenblatt solutions.
A self-similar, radially symmetric and non-negative solution on Ω := R n is called a Barenblatt solution. A Barenblatt solution exists if and only if p(m−1)+n(p−m) > 0, see [17, (12.6)]. In this case it is the solution of the doubly nonlinear diffusion equation (3) to a Dirac measure M δ 0 with a specific mass M > 0 at the origin as in the case p = m [17, (12.9)] resp. by in the case p = m [17, (12.10)] with a constant C > 0 determined by the initial mass M . For m < p Barenblatt solutions have compact support, while for m ≥ p they are positive on the whole space R n .

Derivation of Barenblatt solutions.
To prepare our study of solutions of Barenblatttype for the doubly nonlinear reaction-diffusion equation (1), let us recall how Barenblatt solutions u can be derived for (3), i.e. in the case of pure diffusion. A self-similar, radially symmetric and non-negative function u = 0 has the form u(t, x) = 1 t β y |x| t α with a function y ≥ 0, y = 0. Thus, due to for a Barenblatt solution u the doubly nonlinear diffusion equation (3) reduces to for y in dependence on the variable ξ = t −α |x|. If additionally β(m − 1) = nα, then a further multiplication of (6) by ξ n−1 gives Hence, under the assumption that the terms in the brackets on both sides converge to zero for ξ → ∞, the function y has to solve (y ) p−1 = −αξy m−1 , 3. Solutions of Barenblatt-type in presence of reactions. If the reaction term λu q−1 is present on the right hand side as in (1), then a self-similar and radially symmetric solution In the case m = q, this equation becomes an ODE for y in dependence on the variable In fact, with these α, β the former equation can be multiplied by In the semilinear case m = 2 = p equation (7) is identical with [8, (2.3)] and [12, (2.9)]. A further multiplication of (7) by ξ n−1 gives Remark 1. In the case of pure diffusion only the condition (β +α)(p−1)+α = β(m−1)+1 was necessary to obtain an ODE for y in dependence on the variable ξ. The condition β(m − 1) = nα was merely required for integrability of this ODE. Without this condition (6) is identical with (7) in the case λ = 0, and there is a one-parameter family of (α, β) satisfying the first condition. Particularly, also in the pure diffusive case there exist signchanging self-similar and radially symmetric solutions, as the proof given below is also valid in this case.
In contrast to the case of pure diffusion (8) can not be integrated directly. Moreover, because the coefficient of (y ) p−1 in (7) is singular at ξ = 0, existence of solutions to initial values y(0) = y 0 and y (0) = y 1 can only be expected for y 1 = 0. However, even for such initial values existence of solutions does not follow from the usual theory for ODEs.
To prove local existence of solutions of (7) for initial values y(0) = y 0 and y (0) = 0, let us consider similarly as in [8, Proof of Proposition 3.1] the integral equation corresponding to (8). Due to the boundary conditions y(0) = y 0 and y (0) = 0 equation (8) is equivalent to the integral equation and a further integration gives Because of 0 ≤ s ≤ r we have s n−1 r n−1 ≤ 1, thus for ξ ∈ [0, ξ 0 ] with ξ 0 > 0 sufficiently small the right hand side of the former equation defines a compact operator y → Ky on C([0, ξ 0 ]) with K < 1. Hence, the Leray-Schauder degree deg(Id −K, B 1 (0), 0) of Id −K on the unit ball B 1 (0) ⊂ C([0, ξ 0 ]) is well-defined and identical with deg(Id, B 1 (0), 0) = 1, as the homotopy H(t, y) := y − tKy, t ∈ [0, 1], is admissible due to K < 1, i.e. there exists a solution y ∈ B 1 (0) of y − Ky = 0 or equivalently a local solution of (7) to initial values y(0) = y 0 and y (0) = 0. Note that the former argumentation is called Leray-Schauder fixed point theorem or Rothe's version of the Schauder fixed point theorem [4, Example 5.2.14], and that its application to (1) has already been suggested in [11]. However, to obtain a self-similar solution we do not only need local existence, but global existence of y on [0, ∞). Therefore, note that multiplying (7) by y implies Particularly, the energy decreases in the case α ≥ 0 (where the right hand side is non-positive for ξ > 0), and as boundedness of the energy from above implies boundedness of y and y in the case β > 0, existence on [0, ∞) follows for q ≥ max(m, p). Further, in this case |y(ξ)| ≤ |y 0 | for every ξ > 0 due to y (0) = 0 and E(ξ) ≤ E(0). Similarly, in the case m < q < p (where α < 0, β > 0) and p > 2 the former equation implies the differential inequality p |y | p for p > 2, thus the energy E stays finite for finite ξ < ∞, and again existence on [0, ∞) follows. Therefore, we have proved the following Lemma.
Remark 2. Note that Lemma 3.1 does merely guarantee that y and ξ n−1 (y ) p−1 +αξ n y m−1 are C 1 -functions on (0, ∞), but it does not guarantee that y m−1 is differentiable at points ξ > 0 with y(ξ) = 0 in the case m < 2. Particularly, in this case (7) or the equivalent system may have no sense at points ξ > 0 with y(ξ) = 0. In contrast, if m ≥ 2, then not only y m−1 but due to the validity of (9) also (y ) p−1 is a C 1 -function on (0, ∞), and (7) as well as (12) hold.
Next we claim that for y 0 > 0 the solution y of (8) becomes negative after a finite ξ − > 0, i.e. y is sign-changing.
Proof. Let us consider the non-autonomous planar system associated to (8) via the substitution z := ξ n−1 (y ) p−1 + αξ n y m−1 , which reads as Starting from y(0) = y 0 > 0 and z(0) = 0, it directly follows from the second equation of (13) that z is stricly monotone decreasing as long as y > 0. Further, the first equation of (13) can be rewritten as for y ≥ 0 and z ≤ 0, i.e. y is strictly decreasing as long as y > 0 and z < 0.
The direct analysis of (13) in the first paragraph of the proof of Lemma 3.2 also implies the following corollary.
Let us now focus on the case m < p < q < q c , where α, β > 0. After the solution has entered quadrant I or III, the component y of a solution changes its monotonicity exactly once, namely when the solution curve (y(ξ), z(ξ)) intersects the family of curves given by the zero set of the equation for y in (8). Let us call such a point in quadrant I or III a turning point of the solution, and note that the the zero set of y in (8) is given by the equation i.e. z is just a root (m < 2) or power (m > 2) of y, whose prefactor increases as ξ increases.
Having reached its turning point in quadrant I or III, a solution stays between the ξdependent family of curves given by the zero set (14) of y and the ξ-dependent family of level curves E(ξ, y, z) = E 0 , where E 0 denotes the level of the energy at the turning point. In fact, on the one hand the ξ-dependent vectorfield given by the right hand side of (8) points downwards along the curve (14) in the first quadrant, as y = 0 and z < 0 (resp. the vectorfield points upwards in the third quadrant, as y = 0 and z > 0), thus the solution stays on the right side (resp. left side) of this curve. On the other hand, with the energy E from (11) given in coordinates (ξ, y, z) by and the level E 0 := E(ξ * , y(ξ * ), z(ξ * )) of energy at the turning point (y(ξ * ), z(ξ * )), a solution in the first quadrant (resp. third quadrant) stays after its turning point on the left side (resp. right side) of the ξ-dependent family of curves given implicitly by E(ξ, y, z) = E 0 , because E decreases along solutions, i.e. E is a ξ-dependent Lyapunov function. Hence, after reaching its turning point in quadrant I and III the solution stays in the region between the curves given by (14) and E(ξ, y, z) = E 0 , which is indicated in figure 1. Now from the shape of the zero set (14) of y and of the level set of the energy (15) at the level E 0 it can be seen that for finite ξ there is a gap between the intersection point (0, 0) of (14) with the y-axes and the intersection point (0, y E (ξ)) (y E (ξ) > 0 in case of the first quadrant resp. y E (ξ) < 0 in case of the third quadrant) of the energy level set E(ξ, y, z) = E 0 with the y-axes. Thus, taking additionally into account that z decreases, a solution curve has only two possibilities: Either it escapes the quadrant via the interval [0, y E (ξ)] resp. [y E (ξ), 0] on the y-axes to the fourth resp. second quadrant, or it stays for all ξ > ξ * inside the region between the two families of curves.
However, in both cases the solution converges to the origin as ξ ∞, because due to (10) the ξ-dependent Lyapunov function E is not only decreasing along a solution but strictly decreasing. In fact, the right hand side of (10) is strictly less than zero for y = 0, and y = 0 holds exactly on the family of curves given by (14), but points on this family of curves = (0, 0) are immediately left by every solution (as the ξ-dependent vector field on this curve points downwards resp. upwards). Thus, although the total derivative d dξ E(ξ, y(ξ), z(ξ)) may vanish for certain ξ, the function E(ξ, y(ξ), z(ξ)) is strictly decreasing. Hence, as E attains for every ξ > 0 its minimum precisely at (y, z) = (0, 0), the origin is asymptotically stable. Let us summarize this observation for the ODE (8).
The former lemma gives rise to the question, in which way y converges to 0 as ξ ∞. This seems to be an open problem, and unfortunately we are not able to give a final answer. In fact, although numerics indicate that there are no solutions of (13) to initial values y(0) = y 0 > 0, z(0) = 0, which turn infinitely many times around the origin when converging to (0, 0) as ξ ∞, our former argumentation does not exclude such solutions. Thus, it even seems to be an open problem whether y oscillates as ξ ∞ or not. Let us assume (as numerics indicate) that y does not oscillate but approaches 0 from one side as ξ ∞. Then it would be interesting to know whether y becomes extinct in finite time (note that for pure diffusion the Barenblatt solutions (4) have compact support in the case m < p). Numerically, it seems so that most solutions of (13) tend to (0, 0) along the family of curves (14). Although they never reach this family of curves, let us calculate on the curve z = αξ n y m−1 the solution of the second equation of (13) (as y = 0 on this curve). Choosing y a dependent variable this equation reads as while for z as dependent variable we obtain In both cases the first term on the right hand side determines the behaviour when approaching 0 as ξ ∞. As the equation h = − 1 ξ h has solutions h(ξ) = C 1 ξ , both cases imply convergence to 0 like 1 ξ as ξ ∞. This seems to indicate that there is no extinction in finite time. However, recall that we merely considered the artificial situation, where (y, z) tend to (0, 0) on the family of curves (14) and the second equation of (13) is satisfied, while solutions of (13) never intersect this family of curves again after a turning point in the first or third quadrant. Moreover, solutions to certain initial values even do not seem to approach (0, 0) along the family of curves (14), but along the family of level curves E(ξ, y, z) = E 0 . Thus, the former argumentation does not exclude extinction in finite time.
If y does not become extinct in finite time, then it would be interesting to know the rate of decay of y when approaching 0. In the former artificial situation both y and z converge to 0 like 1 ξ as ξ ∞, but let us again remark that this artificial situation does not need to have anything to do with the real situation. Unfortunately, our argumentation does not allow us to obtain results about the decay rate of y. However, this rate is important when discussing initial values: Let us assume that y either becomes extinct in finite time or tends so fast to 0 that ∞ 0 |y(ξ)| µ−1 ξ n−1 dξ < +∞ holds with µ < m determined by (µ−1)β = nα, e.g. y(ξ) = O( 1 ξ k ) as ξ ∞ with k > p q−p . Then for every test function φ ∈ C c (R n ), i.e. the renormalized initial value u µ−1 0 of the Barenblatttype solution u is a multiple of the Dirac measure δ 0 at the origin.