Recent progress in the theory of Nonlinear Diffusion with Fractional Laplacian Operators

We report on recent progress in the study of nonlinear diffusion equations involving nonlocal, long-range diffusion effects. Our main concern is the so-called fractional porous medium equation, $\partial_t u +(-\Delta)^{s}(u^m)=0$, and some of its generalizations. Contrary to usual porous medium flows, the fractional version has infinite speed of propagation for all exponents $00$; on the other hand, it also generates an $L^1$-contraction semigroup which depends continuously on the exponent of fractional differentiation and the exponent of the nonlinearity. After establishing the general existence and uniqueness theory, the main properties are described: positivity, regularity, continuous dependence, a priori estimates, Schwarz symmetrization, among others. Self-similar solutions are constructed (fractional Barenblatt solutions) and they explain the asymptotic behaviour of a large class of solutions. In the fast diffusion range we study extinction in finite time and we find suitable special solutions. We discuss KPP type propagation. We also examine some related equations that extend the model and briefly comment on current work.


Introduction
This is a follow-up to a previous survey by the author [92], which reported on recent research on two models of nonlinear diffusion processes involving long-range diffusion, written on the occasion of the Abel Symposium held in Oslo in 2010. Such evolution processes are represented by nonlinear parabolic equations involving nonlocal operators of the fractional Laplacian type.
Recapitulating what was said there, the classical theory of diffusion is expressed mathematically by means of the heat equation, and more generally by parabolic equations of linear type; it has had an enormous success and is now a foundation stone in science and technology. The last half of the past century has witnessed intense activity and progress in the theories of nonlinear diffusion, examples being the Stefan Problem, the Porous Medium Equation, the p-Laplacian equation, the Total Variation Flow, evolution problems of Hele-Shaw type, the Keller-Segel chemotaxis system, and many others. Reaction diffusion has also attracted considerable attention. In the last decade there has been a surge of activity focused on the use of so-called fractional diffusion operators to replace the standard Laplace operator (and other kinds of elliptic operators with variable coefficients), with the aim of further extending the theory by taking into account the presence of the long range interactions that occur in a number applications. The new operators do not act by pointwise differentiation but by a global integration with respect to a singular kernel; in that way the nonlocal character of the process is expressed.

Fractional operators
Though there is a wide class of interesting nonlocal operators under scrutiny, a substantial part of the current work deals diffusion modeled by the so-called fractional Laplacians. We recall that the fractional Laplacian operator is a kind of isotropic differentiation operator of order 2s, for some s ∈ (0, 1), that can be conveniently defined through its Fourier Transform symbol, which is |ξ| 2s . Thus, if g is a function in the Schwartz class in R N , N ≥ 1, we write (−∆) s g = h if (1) h (ξ) = |ξ| 2s g(ξ) , so that for s = 1 we recover the standard Laplacian. This definition allows for a wider range of parameters s. The interval of interest for fractional diffusion is 0 < s ≤ 1, and for s < 1 we can also use the integral representation where P.V. stands for principal value and C N,σ is a normalization constant, with precise value C N,2s = 2 2s sΓ((N +2s)/2)/(π N/2 Γ(1−s)). In the limits s → 0 and s → 1 it is possible to recover respectively the identity or the standard minus Laplacian, −∆, cf. [30,72]. Remarkably, the latter one cannot be represented by a nonlocal formula of the type (2). It is also useful to recall that the operators (−∆) −s , 0 < s < 1, inverse of the former ones, are given by standard convolution expressions: in terms of the usual Riesz potentials. Basic references for these operators are the books by Landkof [69] and Stein [83]. A word of caution: in the literature we often find the notation σ = 2s, and then the desired interval is 0 < σ < 2. According to that practice, we will sometimes use σ instead of s.
Another option is to use the classical way of defining the fractional powers of a linear selfadjoint nonnegative operator, and it is expressed in terms of the semigroup associated to such an operator. In the case of the standard Laplacian operator, it reads (4) (−∆) s g(x) = 1 Γ(−s) ∞ 0 e t∆ g(x) − g(x) dt t 1+s .
All the above definitions are equivalent when dealing with the Laplacian on the whole space R N .
The interest in these fractional operators has a long history in Probability for reasons we explain in the next paragraph. Motivation from Mechanics appears in the famous Signorini problem (with α = 1/2), cf. [78,34]. And there are applications in Fluid Mechanics, cf. [41,67]. There is a wide literature on the subject, both for its relevance to Analysis, PDEs, Potential Theory, Stochastic Processes, and Finance, and for the growing number of practical applications. See e.g. [37,54,88,92] where further references can be found.
The systematic study of the corresponding PDE models with fractional operators is relatively recent, and many of the results have been established in the last decade. A part of the current research concerns linear or quasilinear equations of elliptic type. This is a huge subject with well-known classical references that will not be discussed here.

Linear evolution processes
A significant part of the motivation and also a large part of the recent literature is related to evolution problems, that we will discuss in the sequel. Thus, the difference between the standard and the fractional Laplacian is best seen in the stochastic point of view (the stochastic theory of diffusion), and it consists in re-examining the conditions under which the Brownian motion is derived, and taking into account long-range interactions instead of the usual interaction driven by close neighbors. This change of model explains characteristic new features of great importance, like enhanced propagation with the appearance of fat tails at long distances (such tails are to be compared with the typical exponentially small tails of the standard diffusion, or the compactly supported solutions of porous medium flows). Moreover, the space scale of the propagation of the distribution is not proportional to t 1/2 as in the Brownian motion, but to another power of time, that can be adjusted in the model; this is known as anomalous diffusion. The fractional Laplacian operators of the form (−∆) σ/2 , σ ∈ (0, 2), are actually the infinitesimal generators of stable Lévy processes [6,16,44].
We will be concerned with evolution partial differential equations that combine diffusion and fractional operators. As already mentioned, a great variety of diffusive problems in nature, namely those referred to as normal diffusion, are satisfactorily described by the classical Heat Equation or Fokker-Planck linear equation. However, anomalous diffusion is nowadays intensively studied, both theoretically and experimentally since it conveniently explains a number of phenomena in several areas of physics, finance, biology, ecology, geophysics, and many others, which can be briefly summarized as having non-Brownian scaling. This leads to linear anomalous diffusion equations, see e.g. [6,48,100]. Fractional kinetic equations of the diffusion, diffusion-advection, and Fokker-Planck type represent a useful approach for the description of transport dynamics in complex systems which are governed by anomalous diffusion. These fractional equations are usually derived asymptotically from basic random walk models, cf. [61,63,74,73,88,97,99] and their references.
The standard linear evolution equation involving fractional diffusion is a usual model for anomalous diffusion. The equation is solved with the aid of well-known functional analysis tools, like Fourier transform; for instance, it is proved that it generates a semigroup of ordered contractions in L 1 (R N ). Moreover, in this setting it has the integral representation where K s has Fourier transform K s (ξ, t) = e −|ξ| 2s t . This means that, for 0 < s < 1, the kernel K s has the form for some profile function F = F s that is positive and decreasing, and it behaves at infinity like F (r) ∼ r −(N +2s) , [22]. When s = 1/2, F is explicit: If s = 1 the function K s=1 is the Gaussian heat kernel, which has a negative square exponential tail, i.e., a completely different asymptotic behavior.
An integral representation of the evolution of the form (6) is not available in the nonlinear models coming from the applications. This implies the need for new methods, thus motivating our work to be described below.

Nonlinear diffusion models
A main feature of current research in the area of PDEs is the interest in nonlinear equations and systems. There are a number of models of evolution equations that can be considered nonlinear counterparts of the linear fractional heat equation, and combine Laplace operators and nonlinearities in different ways. Let us mention some of the most popular in the recent PDE literature.
• Type I. A natural option is to consider the equation (9) ∂ t u + (−∆) s (u m ) = 0 with 0 < s < 1 and m > 0. This is mathematically the fractional version of the standard Porous Medium Equation (PME) that is recovered as the limit s → 1 and has been extensively studied, cf. [7,91]. We will call equation (9) the Fractional Porous Medium Equation, FPME, as proposed in our first works on the subject [51,52], in collaboration with de A. de Pablo, F. Quirós and A. Rodríguez.
Interest in studying the nonlinear model we propose is two-fold: on the one hand, experts in the mathematics of diffusion want to understand the combination of fractional operators with porous medium type propagation (which is described as degenerate parabolic). Models of this kind arise in statistical mechanics when modeling for instance heat conduction with anomalous properties. On the other hand, it is mentioned in heat control problems by [10]. The rigorous study of such nonlinear models has been delayed until this time by the mathematical difficulties in treating at the same time the nonlinearity and fractional diffusion.
We will devote most of this paper to explain the main results obtained so far on the mathematical theory of this equation. Thus, Section 3 contains the theory of existence, uniqueness and continuous dependence of the Cauchy problem posed in the whole space R N , developed in [51,52]. A main result is the property of infinite speed of propagation. Let us explain the result at this moment and point out the contrast to the PME. While one of the best known features of the PME for m > 1 is the fact that compactly supported nonnegative initial data give rise to solutions that are also compactly supported as functions of x for every fixed t > 0, this property does not hold in equation (9) if s < 1. Indeed, the unique nonnegative solutions we have constructed for nontrivial data u(x, 0) ∈ L 1 + (R N ) are positive everywhere for all t > 0. This happens for all m > 0 and all 0 < s < 1 (for m close to 0 attention must be paid to extinction of the solution in finite time, but that is another issue, see below). Some traits are common with the PME equation: an L 1 -contraction semigroup is constructed and it depends continuously on the exponent of fractional derivation and the exponent of the nonlinearity.
The further investigation of the FPME has taken a number of directions. Thus, in Section 4 we discuss the existence of self-similar solutions with conserved finite mass (we call them fractional Barenblatt solutions), following the study made in [93]. The main properties are derived, allowing to construct the self-similar profiles. Generally, such profiles are not explicit. These solutions are then used to explain the asymptotic behavior of the whole class of nonnegative solutions with finite mass (a kind of fractional central limit theorem).
The study of a priori estimates (some of them universal), quantitative bounds on positivity and Harnack estimates is done in Section 5 following the results obtained in [27] in collaboration with M. Bonforte. This is done in the favorite setup, the Cauchy problem. But the fractional Laplace operators offer many novelties when posed in a bounded domain of the Euclidean space, say with Dirichlet boundary conditions. We report in Section 6 on recent work done with M. Bonforte on the issue and comment on related works.
As a natural extension of the equation, we have recently considered equations of the form where Φ is a monotone increasing function. Collaboration with Pablo, Quirós and Rodríguez, [53] and [94]; regularity is a main issue, we review the results in Section 7. Section 8 reports on the technique of Schwarz symmetrization of the solutions of the these problems and derives a priori estimates. This is the result of collaboration with B. Volzone [95,96]. Finally, numerical methods for these equations are being investigated by our team and a number of authors, and we will make a comment on that issue.
• The combination of fractional diffusion with other effects leads to interesting behaviour. One of the most common combinations is reaction-diffusion. In that sense there has been interesting work on the possible extension of the well-known KPP behaviour to cover fractional diffusion, both in the linear and nonlinear case. Results by X. Cabré and J. M. Roquejoffre [33] (for linear diffusion) and Diana Stan and the author [81] (in the nonlinear case) are covered in Section 9.
On the other hand, the combination of fractional diffusion convection appears in a number of models of the so-called geostrophic equations, like (12) ∂ t u + v · ∇u = (−∆) s u, where 0 < s < 1 and v(x, t) is a divergence free vector field related to u in different ways. Important works in this direction related to our research are due to Caffarelli and Vasseur [41] and Kiselev et al. [67], but again, we will not deal with this topic here. This last model uses linear fractional diffusion. Nonlinear diffusion combined with convection has been studied by some authors, like Cifani and Jakobsen [45], where references to other models are given.
• Type II. We recall here the other model of diffusion with fractional operators where the author has been involved, mainly in collaboration with Luis Caffarelli. Since it was extensively reported in the survey paper [92] we will only give a short reminder. This alternative model is derived in a more classical way from the Porous Medium Equation since it is based on the usual Darcy law, with the novelty that the pressure is related to the density by an inverse fractional Laplacian operator. This nonlinear fractional diffusion equation of porous medium type takes the form where K is the Riesz operator that typically expresses the inverse to the fractional Laplacian, Ku = (−∆) −s u. This has been studied by Caffarelli and the author in [39,40] and Biler, Karch and Monneau [21,19], where the equation is derived in the framework of the theory of dislocations. We recall that this model has some strikingly different properties, like lack of strict positivity and occurrence of free boundaries. Self-similar solutions also exist but their existence and properties are quite different from those of equation (9). Even the asymptotic behaviour is quite different. In order to distinguish both models the name porous medium equation with fractional pressure has been proposed for equation (13).
As for recent work, the boundedness and Hölder regularity of finite mass solutions is studied in joint work with L. Caffarelli and F. Soria, [38]. As a limit case of this second model s → 1, one obtains a variant of the equation for the evolution of vortices in superconductivity derived heuristically by Chapman-Rubinstein-Schatzman [47] and W. E [56] as the hydrodynamic limit of Ginzburg Landau, and studied by Lin and Zhang [71], and Ambrosio and Serfaty [4]. The understanding of this limit has been done in collaboration with Sylvia Serfaty [76], and is related to work by Bertozzi et al. on aggregation models [17,18].
The more general equation u t = ∇(u ∇K(u m−1 )), m > 0, is considered in [20], and self-similar solutions are found in explicit form, a very interesting fact. A very natural extension of these equations is (14) u where we may take for simplicity u ≥ 0, m > 1 and p > 0. We are working in such generalizations, cf. [82] where the property of finite speed of propagation is examined.
The study of the fine asymptotic behavior, i.e., obtaining rates of convergence is not easy. Progress on this issue is under way in collaboration with Carrillo and Huang [43].
• Disclaimer. Since our objective is limited, we are leaving completely outside of the presentation many related topics. Thus, the anomalous diffusion that is encountered in many applications is often described by means of fractional time operators, using the so-called Riemann-Liouville derivatives. There is a huge literature on the subject in the case of linear equations. We will not enter into such topic in this document, though some nonlinear models are quite appealing.

The fractional porous medium equation
We now turn our attention to the nonlinear heat equation with fractional diffusion (9). Indeed, it is a whole family of equations with exponents s ∈ (0, 1) and m > 0. They can be seen as fractional-diffusion versions of the PME described above, [91], [90]. The classical Heat Equation is recovered in this model in the limit s = 1 when m = 1, the PME when m > 1, the Fast Diffusion Equation when m < 1.

Some applied literature
We gather here some updated information on the occurrence of the nonlinear fractional diffusion equation we propose and related models in the physical or probabilistic literature.
• Anomalous diffusion often takes a nonlinear form. To be more specific, there exist many phenomena in nature where, as time goes on, a crossover is observed between different diffusion regimes. Tsallis et al. [23,70] discuss the following cases: (i) a mixture of the porous medium equation, which is connected with non-extensive statistical mechanics, with the normal diffusion equation; (ii) a mixture of the fractional time derivative and normal diffusion equations; (iii) a mixture of the fractional space derivative, which is related with Lévy flights, and normal diffusion equations. In all three cases a crossover is obtained between anomalous and normal diffusions. This leads to models of nonlinear diffusion of porous medium or fast diffusion types with standard or fractional Laplace operators, cf. equation (4) of [23].
• There have been many studies of hydrodynamic limits of interacting particle systems with long-range dynamics, which lead to fractional diffusion equations of our type, mainly linear like in [63], but also nonlinear in the recent literature, cf. the works [62], [64]. Thus, in the last reference, Jara and co-authors study the non-equilibrium functional central limit theorem for the position of a tagged particle in a mean-zero one-dimensional zero-range process. The asymptotic behavior of the particle is described by a stochastic differential equation governed by the solution of the following nonlinear hydrodynamic (PDE) equation, ∂ t ρ = a 2 ∂ 2 x Φ(ρ). When Φ is a power we recover equation (9).
• Equations like the last one (in several space dimensions) occur in boundary heat control, as already mentioned by Athanasopoulos and Caffarelli [11] , where they refer to the model formulated in the book by Duvaut and Lions [55], and use the so-called Caffarelli-Silvestre extension [37].

Mathematical problem and general notions
Let us present the main features and results in the theory we have developed. To be specific, the theory of existence and uniqueness as well the main properties are first studied by De Pablo, Quirós, Rodríguez, and Vázquez in [51,52] for the Cauchy problem We have put σ = 2s. The notation |u| m−1 u is used to allow for solutions of two signs, but since it is a bit awkward we will often write u m instead of |u| m−1 u even for signed solutions when no confusion is feared; the same happens with the power u 1/m . We take initial data f ∈ L 1 (R N ), which is a standard assumption in diffusion problems on physical grounds. As for the exponents, we consider the whole fractional exponent range 0 < σ < 2, and take porous medium exponent m > 0. As we have said, in the limit σ → 2 we want to recover the standard The two papers contain a rather complete analysis of the problem. A semigroup of weak energy solutions is constructed for every choice of m and σ, both the L ∞ smoothing effect and C α regularity work in most cases (with some restrictions if m is near 0), and there is infinite propagation for all m > 0 and σ < 1. The results can be viewed as a nonlinear interpolation between the extreme cases σ = 2: u t − ∆(|u| m−1 u) = 0, and σ = 0 which turns out to be a simple ODE: u t + |u| m−1 u = 0. It is to be noted that the critical exponent m c := (N − σ) + /N plays a role in the qualitative theory: the properties of the semigroup are more familiar when m > m c . A similar exponent is well-known in the Fast Diffusion theory (putting σ = 2). Note that such exponent is not considered when N = 1 and σ ≥ 1.
Preliminary notions. If ψ and ϕ belong to the Schwartz class, the definition (1) of the fractional Laplacian together with Plancherel's theorem yield Therefore, if we multiply the equation in (15) by a test function ϕ and integrate by parts, we obtain This identity is the basis of our definition of a weak solution. The integrals in (16) make sense if u and u m belong to suitable spaces. The right space for u m is the fractional Sobolev spacė H σ/2 (R N ), defined as the completion of C ∞ 0 (R N ) with the norm where ψ denotes the Fourier transform of ψ. Note that Definition A function u is a weak solution to Problem (15) if: A drawback of this definition is that there is no convenient formula for the fractional Laplacian of a product or of a composition of functions. Moreover, we take no advantage in using compactly supported test functions since their fractional Laplacian loses this property. To overcome these and other difficulties, we will use the fact that our solution u is the trace of the solution of a local problem obtained by extending u m to a half-space whose boundary is our original space.
Extension Method. In the particular case σ = 1 studied in [51], the problem is reformulated by means of the well-known representation of the half-Laplacian in terms of the Dirichlet-Neumann operator. This allowed us to transform the nonlocal problem into a local one (i. e., involving only derivatives and not integral operators). Of course, this simplification pays a prize, namely, introducing an extra space variable. The application of such an idea is not so simple when σ = 1; it involves a number of difficulties that we address in [52]. We have to use the characterization of the Laplacian of order σ, (−∆) σ/2 , 0 < σ < 2, described by Caffarelli and Silvestre in their famous paper [37], in terms of the so-called σ-harmonic extension, which is the solution of an elliptic problem with a degenerate or singular weight.
Let us explain this extension in some more detail. Then, where the precise constant, which does not depend on N , is µ σ = 2 σ−1 Γ(σ/2) Γ(1−σ/2) , see [37]. In (17) the operator ∇ acts in all (x, y) variables, while in (18) (−∆) σ/2 acts only on the x = (x 1 , · · · , x N ) variables. In the sequel we denote Notation. In this section we will use the notation Ω = R N × (0, ∞) for the upper half-space, with points x = (x, y), x ∈ R N , y > 0; its boundary, which is identified to the original R N with variable x, will be named Γ. Besides, we use the simplified notation u m for data of any sign, instead of the actual "odd power" |u| m−1 u, and we will also use such a notation when m is replaced by 1/m. The convention is not applied to any other powers.
Extended problem. Weak solutions. With the above in mind, we rewrite problem (15) for w = u m as a quasi-stationary problem with a dynamical boundary condition This problem has been considered by Athanasopoulos and Caffarelli [11], who prove that any bounded weak solution is Hölder continuous if m > 1.
To define a weak solution of this problem we multiply formally the equation in (19) by a test function ϕ and integrate by parts to obtain where u = (Tr(w)) 1/m is the trace of w on Γ to the power 1/m. This holds on the condition that ϕ vanishes for t = 0 and t = T , and also for large |x| and y. We then introduce the energy space X σ (Ω), the completion of C ∞ 0 (Ω) with the norm The trace operator is well defined in this space, see below.
Definition A pair of functions (u, w) is a weak solution to Problem (19) if: For brevity we will refer sometimes to the solution as only u, or even only w, when no confusion arises, since it is clear how to complete the pair from one of the components, u = (Tr(w)) 1/m , w = E(u m ).
Equivalence of weak formulations. The key point of the above discussion is that the definitions of weak solution for our original nonlocal problem and for the extended local problem are equivalent. The main ingredient of the proof is that equation (18) holds in the sense of distributions for any g ∈Ḣ σ/2 (Γ).
Proposition A function u is a weak solution to Problem (15) if and only if (u, E(u m )) is a weak solution to Problem (19).
Strong solutions. Weak solutions satisfy equation (15) in the sense of distributions. Hence, if the left hand side is a function, the right hand side is also a function and the equation holds almost everywhere. This fact allows to prove uniqueness and several other important properties, and hence motivates the following definition.

Main results
Existence. We prove existence of a suitable concept of (weak) solution for general L 1 initial data only in the restricted range m > m c , which includes as a particular case the linear fractional heat equation, case m = 1. If 0 < m ≤ m c (which implies that 0 < σ < 1 if N = 1, recall that m c = (N − σ) + /N ) we need to slightly restrict the data to obtain weak solutions.
there exists a weak solution to the Cauchy problem for the FPME.
Uniqueness. We first prove uniqueness of solutions in the range m ≥ m c . If 0 < m < m c , we need to use the concept of strong solution, a concept that is standard in the abstract theory of evolution equations. This is no restriction in view of the next results proved in [52]. We state the uniqueness result in its simplest version. Qualitative behaviour. The solutions to Problem (15) have some nice properties that are summarized here.
Theorem 3.4. Assume f, f 1 , f 2 satisfy the hypotheses of Theorem 3.1, and let u, u 1 , u 2 be the corresponding strong solutions to Problem (15).
In the linear case m = 1 the above properties: conservation of mass, the smoothing effect with a precise decay rate, positivity and regularity, can be derived directly from the representation formula (6) and the properties of the kernel K σ .
Continuous dependence. We also show that the solution (i.e., the semigroup) depends continuously on the initial data and on both parameters m and σ, in particular in the nontrivial limit σ → 2, that allows to recover the standard PME, ∂ t u − ∆|u| m−1 u = 0, or the other end σ → 0, for which we get the ODE : ∂ t u + |u| m−1 u = 0. Continuity will be true in general only in L 1 loc , unless we stay in the region of parameters where mass is conserved. Theorem 3.5. The strong solutions depend continuously in the norm of the space C([0, T ] : L 1 loc (R N )) on the parameters m, σ, and the initial data f . If moreover m ≥ m c and 0 < σ ≤ 2, convergence also holds in C([0, T ] : L 1 (R N )).
We will extend this theory in different ways. See in particular Section 7.

Self-similar solutions and their role
Several types of special solutions play a role in the theory of the FPME: fundamental solutions, very singular solutions, extinction solutions,...

Fundamental solutions
Let us consider the Cauchy problem for the FPME taking as initial data a Dirac delta, (23) u Solutions with such data are called fundamental solutions in the linear theory, and we will keep that name though their relevance is different in the nonlinear context. We use the concept of continuous and nonnegative weak solution introduced in [51], [52], for which there is a welldeveloped theory when the datum is a function in L 1 (R N ) as we have outlined. We recall that the question under discussion has a well-known answer for the standard Porous Medium Equation, i.e., the limit case s = 1, in the form of the usual Barenblatt solutions; they were actually discovered in the 1950's by Zeldovich-Kompanyeets [101] and Barenblatt [14], cf. their use in applications in [15] and in theory of the equation in [91].
Let us describe the results of our paper [93], devoted to prove existence, uniqueness and main properties of the fundamental solutions for the FPME. The exponent m varies in principle in the range m > 1, but the methods extend to the linear case m = 1 and even to the fast diffusion range (FDE) m < 1 on the condition that m > m c = max{(N − 2s)/N, 0}. Such a type of restriction on m from below carries over from the PME-FDE theory, [90]. The solutions we find are self-similar functions of the form with suitable exponents α and β and profile F ≥ 0. Here is the main result.
The profile function F M (r), r ≥ 0, is a bounded and continuous function, it is positive everywhere, it is monotone and it goes to zero at infinity.
Let us point out that the Brownian scaling (β = 1/2 and α = N/2) is recovered in the fractional case under the condition N (m − 1) = 2(1 − s), for instance for m = (N + 1)/N in the most common case s = 1/2. Moreover, the profile F M can be obtained by a simple rescaling of F 1 . F is actually a smooth function by the regularity results of [94].
The initial data are taken in the weak sense of measures (26) lim for all φ ∈ C b (R N ), the space of continuous and bounded functions in R N . We will call these self-similar solutions of Problem (9)-(23) with given M > 0 the Barenblatt solutions of the fractional diffusion model by analogy with the PME and other prominent studied cases. The form of the exponents explains the already mentioned restriction on m from below.
Idea of the proofs: existence of the Barenblatt solutions is done by approximation of the Cauchy problem with smooth initial data and passing to the limit, using suitable a priori estimates. Uniqueness needs to go over to the potential equation as follows: we take the convolution of u(x, t) with the Riesz kernel and define: Then (−∆) s U = u, and using the equation we get equation (PE): in other words, U t + ((−∆) s U ) m = 0. Though awkward-looking, this 'dual equation' admits a good uniqueness theorem, as shown in [93].
The properties of the profile function F (r) are quite important, in particular the asymptotic behavior. This is why in [93] we derive the elliptic equation that it satisfies:   The case m = m 1 is borderline and has a logarithmic correction.
Solutions with explicit form. As this survey is written, Y. Huang reports [60] the explicit expression of the Barenblatt solution for the special value of m, m ex = (N + 2 − 2s)/(N + 2s). The profile is given by where the two constants λ and R are determined by the total mass M of the solution and the parameter β. Note that for s = 1/2 we have m ex = 1, and the solution is the one mentioned in the introduction for the linear case, (8). We always have m ex > m 1 ; the range of m that is covered for 0 < s < 1 is N/(N + 2) < m ex < (N + 2)/N .

Asymptotic behavior
As a main application of this construction, we prove in [93] that the asymptotic behaviour as t → ∞ of the class of nonnegative weak solutions of (9) with finite mass (i. e., R N u(x, t) dx < ∞) is given in first approximation by the family u * M (x, t), i. e., the Barenblatt solutions are the attractors in that class of solutions.
and also uniformly in x ∈ R N . It follows that for every p ∈ (1, ∞) we have This result is a clear example of the important role that the fundamental solutions and their properties can play in the applications.

Very singular solutions
Only in the range m c < m < m 1 = N/(N + 2s) we can pass to the limit M → ∞ and obtain a special solution that has a fixed isolated singularity at x = 0 that has separate-variables form and we call very singular solution (VSS) by analogy with the standard Fast Diffusion Equation. The result represents a marked difference with the case s = 1 where VSS exist in the larger range m c < m < 1. It is another manifestation of the long-range interactions of the fractional Laplacian, that avoids some of the purely local estimates of the standard FDE with the classical Laplacian operator; indeed, extrapolation of the standard estimates would justify the existence of a VSS in cases where there is none.
Among other interesting qualitative properties of the equation, in [93] we prove an Aleksandrov reflection principle. This is related to symmetrization principles that we will discuss below.

Other solutions
Self-similar solutions of the second kind, where the similarity exponents α and β are only related by the compatibility condition (m − 1)α + 2sβ = 1, but are otherwise free not to verify the constant-mass condition α = N β are constructed in the work with Volzone [96].
The question of extinction of mass in finite time that was introduced in [52] and happens for m < m c < 1 is also best understood in terms of explicit solutions, which were constructed in [96] and have the form (35) U (x, t) 1−m = C T − t |x| 2s .

Estimates, better existence and positivity
In collaboration with M. Bonforte we have studied the existence of quantitative a priori estimates of a local type for solutions of the FPME and we have then derived consequences for the theory. Thus, in [27] we have dealt with the properties of nonnegative solutions of the Cauchy problem. Such kind of estimates were obtained for the standard PME by Aronson-Caffarelli [8] and by the present authors for the standard FDE [24,25] . The same local tools are not efficient for the present fractional model due to the nonlocal character of the diffusion operator, but then estimates occur in weighted spaces.

Weighted L 1 estimates in the fast diffusion range
We will concentrate on the case m < 1. The use of suitable weight functions allows to prove crucial L 1 -weighted estimates that enter substantially into the derivation of the main results. The results take different forms according to the value of the exponent m, a fact that is to be expected since it happens for standard FDE (i.e., the limit case s = 1). When s < 1 the equation is nonlocal, therefore we cannot expect purely local estimates to hold. Indeed, we will obtain estimates in weighted spaces if the weight satisfies certain decay conditions at infinity. We present first a technical lemma that shows the difference with the standard Laplacian.
Theorem 5.2. Let u ≥ v be two ordered solutions to the FPME with 0 < m < 1. Let ϕ R (x) = ϕ(x/R) where R > 0 and ϕ is as in the previous lemma with 0 ≤ ϕ(x) ≤ |x| −α for |x| >> 1 and Then, for all 0 ≤ τ, t < ∞ we have It is remarkable that the estimate holds for (very) weak solutions, even changing sign solutions. Also, it is worth pointing out that the estimate holds both for τ < t and for τ > t. The estimate implies the conservation of mass when m c < m < 1, by letting R → ∞. On the other hand, when 0 < m < m c solutions corresponding to u 0 ∈ L 1 (R N ) ∩ L p (R N ) with p ≥ N (1 − m)/2s , extinguish in finite time T (u 0 ) > 0 , (see e.g. [52]); the above estimates provide a lower bound for the extinction time in such a case, just by letting τ = T and t = 0 in the above estimates: Moreover, if the initial datum u 0 is such that the limit as R → +∞ of the right-hand side diverges to +∞, then the corresponding solution u(t, x) exists (and is positive) globally in time.

Existence of solutions in weighted L 1 -spaces
As a first consequence of these estimates, we can extend the L 1 existence theory of [51,52] to an existence result for very weak solutions with non-integrable data in some weighted L 1 ϕ space. In particular, bounded initial data or data with slow growth at infinity are allowed.
. This solution is continuous in the weighted space, u ∈ C([0, T ] : Idea of the proof: We take 0 ≤ u 0,n ∈ L 1 (R N ) ∩ L ∞ (R N ) be a non-decreasing sequence of initial data u 0,n−1 ≤ u 0,n , converging monotonically to u 0 ∈ L 1 (R N , ϕ dx) , i. e., such that R N (u 0 − u n,0 )ϕ dx → 0 as n → ∞. Consider the unique solutions u n (t, x) of the equation with initial data u 0,n . By the comparison results of [52] we know that they form a monotone sequence. The weighted estimates (36) show that the sequence is bounded in L 1 (R N , ϕ dx) uniformly in t ∈ [0, T ] . By the monotone convergence theorem in L 1 (R N , ϕ dx), we know that the solutions u n (t, x) converge monotonically as n → ∞ to a function u(t, x) ∈ L ∞ ((0, T ) : L 1 (R N , ϕ dx)). We then pass to the limit n → ∞.
Remark. The solutions constructed above only need to be integrable with respect to the weight ϕ, which has a tail of order less than d + 2s/m. Therefore, we have proved existence of solutions corresponding to initial data u 0 that can grow at infinity as |x| (2s/m)−ε for any ε > 0 . Note that for the linear case m = 1 this exponent is optimal in view of the representation of solutions in terms of the fundamental solution, but this does not seem to be the case for m < 1.
Theorem 5.4 (Uniqueness). The solution constructed in Theorem 5.3 by approximation from below is unique. We call it the minimal solution. In this class of solutions the standard comparison result holds, and also the estimates of Theorem 5.2.

Lower bounds for fractional fast diffusion
Section 3 of paper [27] studies the actual positivity of nonnegative solutions via quantitative lower estimates in the so called good fast diffusion range. We will use the previous notation: The following is the main quantitative estimate from below for positive solutions.
Let u(t, ·) ∈ L 1 (R N , ϕ dx) be a very weak solution to the FPME with initial datum u 0 . Then can calculate a time (38) t with quantified C * such that The positive constants C * , K 1 , K 2 depend only on m, s and N ≥ 1.
We point out that the diffusion in the equation is nonlocal, but the estimates are local. These estimates can be combined with the L ∞ bounds from above of the papers [51,52] to sandwich the solution from above and below. The paper continues to obtain sharp estimates on the behaviour as |x| → ∞ (so-called tail behaviour).

Estimates in other ranges
New local estimates are obtained in the paper either for 0 < m < m c or for m ≥ 1. The question of Harnack inequalities is also discussed. We refrain from giving more details at this moment for lack of space, and refer to the paper [27] for further information.
An interesting result is the existence and uniqueness of an initial trace for nonnegative solutions that we state in complete detail for m < 1.
Theorem 5.6. Let 0 < m < 1 and let u be a nonnegative weak solution of equation (9) in (0, T ] × R N . Assume that u(T ) L 1 (R N ) < ∞. Then there exists a unique nonnegative Radon measure µ as initial trace, that is Moreover, the initial trace µ satisfies the bound The result is indeed true for m > 0, including of course the linear equation for m = 1.

Problems in bounded domains
The investigation of the questions of existence of solutions for the FPME posed in a bounded domain with suitable boundary conditions was initiated in the papers [51,52]. The study of a priori estimates has been taken up in ongoing collaboration with M. Bonforte. These estimates are quite different from the problem in the whole space in statements and methods.

Fractional Laplacian operators on bounded domains
The first problem that we encounter is the possibility of various reasonable definitions of the concept of fractional Laplacian operator. This is in contrast with the situation in the whole Euclidean space R N , where there is a natural concept of fractional Laplacian that can be defined in several equivalent ways as we have mentioned in Subsection 1.1. When we consider equation (9) posed on bounded domains, the definition via the Fourier transform does not apply, and different choices appear as possible definitions of the fractional Laplacian.
• On one hand, starting from the classical Dirichlet Laplacian ∆ Ω on the domain Ω , the so-called spectral definition of the fractional power of ∆ Ω uses the formula in terms of the semigroup associated to the Laplacian, namely We will call the operator defined in such a way the spectral fractional Laplacian, SFL for short, and denote as L 1 = (−∆ Ω ) s . In this case, the initial and boundary conditions associated to the fractional diffusion equation (9) read (44) u(t, x) = 0 , in (0, ∞) × ∂Ω , u(0, ·) = u 0 , in Ω , Let us list some properties of the operator: L 1 = (−∆ Ω ) s is a self-adjoint operator on L 2 (Ω) , with a discrete spectrum: its eigenvalues are the family of s-power of the eigenvalues of the Dirichlet Laplacian: (λ j ) s > 0, j = 1, 2, . . .; the corresponding normalized eigenfunctions Φ j are exactly the same as the Dirichlet Laplacian, therefore they are as smooth as the boundary allows, namely when ∂Ω is C k , then Φ j ∈ C ∞ (Ω) ∩ C k (Ω) . We can thus write (45) ( The definition of the Fractional Laplacian via the Caffarelli-Silvestre extension [37] has been extended to the case of bounded domains by Cabré and Tan [31] by using as extended domain the cylinder C = (0, ∞) × Ω in R N +1 + , and by putting zero boundary conditions on the lateral boundary of that cylinder. This definition enables to understand the boundary conditions in an easy way. It is proved that this definition is equivalent to the SFL. See also [42].
• On the other hand, we can define a fractional Laplacian operator by using the integral representation (2) in terms of hypersingular kernels and "restrict" the operator to functions that are zero outside Ω: we will denote the operator defined in such a way as L 2 = (−∆ |Ω ) s , and call it the restricted fractional Laplacian, RFL for short. In this case, the initial and boundary conditions associated to the fractional diffusion equation (9) read See more in [28]. We also refer to [77] for a discussion and references about the differences between the Spectral and the Restricted fractional Laplacian. The authors of [77] call the second type simply Fractional Laplacian, but we feel that the absence of descriptive name leads to confusion.

Main results on the Dirichlet problem
In our paper [28] we choose to work with the Dirichlet spectral laplacian, L 1 . For a quite general class of nonnegative weak solutions to the above problem, we derive in absolute upper estimates up to the boundary of the form ∀t > 0 , ∀x ∈ Ω .
In particular, we observe that the boundary behaviour is dictated by Φ 1 , the first positive eigenfunction of L 1 , which behaves like the distance to the boundary at least when the domain is smooth enough. We will also prove standard and weighted instantaneous smoothing effects of the type where ϑ 1,1 = 1/(2s + (N + 1)(m − 1)) . This is sharper than (47) only for small times. As a consequence of the above upper estimates, we derive a number of useful weighted estimates and we also obtain backward in time smoothing effect of the form which is quite surprising and has not been observed before to our knowledge.
We then pass to the question of lower estimates, the main result being the estimate where the waiting time t * has the explicit form . Then we observe that the above estimates combine into global Harnack inequalities in the following form: for all t ≥ t * We also provide as a corollary the local Harnack inequalities of elliptic type.
All the constants in the above results are universal, in the sense that may depend only on N, m, s and Ω , but not on u. They also have an almost explicit expression, usually given in the proof. Actually, we have tried to obtain quantitative versions of the estimates with indication of the dependence of the relevant constants. In some cases the estimates are absolute, in the sense that they are valid independently of the (norm of the) initial data.
New method. It is worth mentioning that the proofs are based on new ideas with respect to what was used in the standard nonlinear diffusion or in the previous section for fractional diffusion in the whole space (both in our papers and in the references). Thus, we exploit the functional properties of the linear operator as much as possible. More precisely, we use the Green function of the fractional operator even in the definition of solution, and we make use of estimates on its behaviour in the proofs. A key ingredient is thus the knowledge of good estimates for the Green function, that we discuss at length in the paper.
A careful inspection of the proofs shows that the presented method would allow to treat a quite wide class of linear operators, an issue that we shall discuss in a forthcoming paper [29]. Here, we have written everything referring to the concrete case of the spectral fractional Laplacian in order to keep the exposition clear and to focus on the main ideas, but the arguments are devised in view of the wider applicability.

Existence theory
The paper is complemented with a brief presentation of the theory of existence and uniqueness of the class of very weak solutions that we use. This complements the basic theory developed before in [51,52] and then extended in [27]. We have no space to enter this intriguing issue here.

The work on general linearities
In collaboration with Pablo, Quirós and Rodríguez we have examined the question of regularity of the nonnegative solutions of the FPME. Since we had proved continuity and positivity, the equation is not more degenerate, at least at the formal level, and the solutions should be smooth if we recall what happens in the classical heat equation and fast diffusion case when s = 1. Actually, in [94] we have decided to study the question for the more general class of equations (52) ∂ t u + (−∆) σ/2 ϕ(u) = 0 .
To be specific, we deal with the Cauchy problem posed in Q = R N × (0, ∞). The constitutive function ϕ is assumed to be at least continuous and nondecreasing. Further conditions will be introduced as needed. This generality was motivated by previous work on the model of logarithmic porous medium type equation with fractional diffusion which is described in detail in [53] and has interesting connections with drift-diffusion models.
The existence of a unique weak solution to the Cauchy problem for the FPME has been fully investigated in [51,52] for the case where ϕ is a positive power. This theory is extended in [94] under suitable conditions on ϕ that cover powers in particular.
Regarding regularity, for the linear fractional heat equation it follows from explicit representation with a kernel that solutions are C ∞ smooth and bounded for every t > 0, x ∈ R N , under the assumption that the initial data are integrable. In the nonlinear case such a representation is not available. Nevertheless, we will still be able to obtain that the solution is smooth if the equation is "uniformly parabolic", 0 < c ≤ ϕ ′ (u) ≤ C < ∞. Our first result establishes that if the nonlinearity is smooth enough, compared to the order of the equation, max{1, σ}, then bounded weak solutions are indeed classical solutions.
The precise regularity of the solution is determined by the regularity of the nonlinearity ϕ; see the paper for the details. Notice that the condition ϕ ′ > 0 together with the boundedness of u implies that the equation is uniformly parabolic.
The idea of the proof is as follows: thanks to the results of Athanasopoulos and Caffarelli [11], we know that bounded weak solutions are C α regular for some α ∈ (0, 1). In order to improve this regularity we write the equation (52) as a fractional linear heat equation with a source term. This term is in principle not very smooth, but it has some good properties. To be precise, given (x 0 , t 0 ) ∈ Q, we have after the time rescaling t → t/ϕ ′ (u(x 0 , t 0 )). A very delicate study of the linear theory, with special attention to the properties of the kernel of the fractional Laplacian, allows to show that solutions to the linear equation (53) have the same regularity as f . Next, using the nonlinearity we observe that f in the actual right-hand side is more regular than u near (x 0 , t 0 ). We are thus in a situation that is somewhat similar to the one considered by Caffarelli and Vasseur in [41], where they deal with an equation, motivated by the study of geostrophic equations, of the form where v is a divergence free vector. Comparing with (53), we see two differences: in their case σ = 1, and the source term is local. These two differences will significantly complicate our analysis.
Singular and degenerate equations. The hypotheses made in Theorem 7.1 excludes all the powers ϕ(u) = |u| m−1 u for m > 0, m = 1, since they are degenerate (m > 1) or singular (m < 1) at the level u = 0. Nevertheless, a close look at our proof shows that we may in fact get a "local" result, see the paper. Therefore, we get for these nonlinearities (and also for more general ones) a regularity result in the positivity (negativity) set of the solution that implies that bounded weak solutions which are either positive or negative are actually classical.
Higher regularity. If ϕ is C ∞ we prove that solutions are C ∞ . The result will be a consequence of the regularity already provided by Theorem 7.1 plus a special result for linear equations with variable coefficients. The case σ < 1 is even more involved since we first have to raise the regularity in space exponent from σ to 1.
Theorem 7.2. Let u be a bounded weak solution to equation (52).
Comments and extensions. Work on extending such results to problems in bounded domains, and to more general operators is in project. Regarding related literature, let us remark that Kiselev et al. [67] give a proof of C ∞ regularity of a class of periodic solutions of geostrophic equations in 2D with C ∞ data. Their methods are completely different to the ones used by us. On the other hand, Cifani and Jakobsen propose in [45] an alternative L 1 theory dealing with a more general class of nonlocal porous medium equations, including strong degeneration and convection. The quantitative study of continuous dependence has been taken up recently in work of Alibaud et al. [1,2].

Estimates via symmetrization
In two papers [95,96] in collaboration with B. Volzone we have used the techniques of Schwarz and Steiner symmetrization to obtain a priori estimates, in many cases with best constants, for the solutions of the FPME, or its more general version u t + (−∆) s A(u) = 0, posed in the whole space or in a bounded domain. The main results concern the case of power nonlinearities, the equation is posed in the whole space and the exponent m < 1. Elliptic results are obtained in [95] as a preliminary to the parabolic theory. In order to keep up with previous sections we will use the letter ϕ for the nonlinearity instead of the A in the paper.
Motivation. Symmetrization is a very old geometrical idea that has become nowadays a popular tool of obtaining a priori estimates for the solutions of different partial differential equations, notably those of elliptic and parabolic type. The application of Schwarz symmetrization to obtaining a priori estimates for elliptic problems is already described in [98].
The standard elliptic result refers to the solutions of an equation of the form posed in a bounded domain Ω ⊆ R N ; the coefficients {a ij } are assumed to be bounded, measurable and satisfy the usual ellipticity condition; finally, we take zero Dirichlet boundary conditions on the boundary ∂Ω. The classical analysis introduced by Talenti [84] leads to pointwise comparison between (i) the symmetrized version (more precisely the spherical decreasing rearrangement) of the actual solution of the problem u(x) and (ii) the radially symmetric solution v(|x|) of some radially symmetric model problem which is posed in a ball with the same volume as Ω. Sharp a priori estimates for the solutions are then derived. Extensions of this method to more general problems or related equations have led to a copious literature.
Parabolic version. This pointwise comparison fails for parabolic problems and the appropriate concept is comparison of concentrations, cf. Bandle [13] and Vázquez [89]. The latter considers the evolution problems of the form where ϕ a monotone increasing real function and u 0 is a suitably given initial datum which is assumed to be integrable. For simplicity the problem was posed for x ∈ R N , but bounded open sets can be used as spatial domains.
Fractional operators. Symmetrization techniques were first applied to PDEs involving fractional Laplacian operators in the paper [24], where the linear elliptic case is studied. In our first paper [95] we were able to improve on that progress and combine it with the parabolic ideas of [89] to establish the relevant comparison theorems based on symmetrization for linear and nonlinear parabolic equations. To be specific, we deal with equations of the form Let us describe the results in [95] concerning the idea of concentration comparison in the case of the solutions to the Cauchy problem Here, the nonlinearity ϕ(u) is a nonnegative function, smooth on R + , with ϕ(0) = 0 and ϕ ′ (u) > 0 for all u > 0 (extended anti-symmetrically in the general two-signed theory). Special attention is paid to cases of the form ϕ(u) = u m with m > 0.In [95] we have obtained that a concentration comparison for solutions to (56) holds only when the nonlinearity ϕ is a concave function, while for convex ϕ a remarkable example is constructed, showing that a failure of concentration comparison occurs (see [95]).
Moreover, the following corollary justifies a reasonable consequence: if the data of problem (56) are less concentrated that those of the symmetrized problem, so are the corresponding solutions.
Corollary 1. With the same assumptions of Theorem 8.1, suppose that u is the solution to problem (56) and v solves where f ∈ L 1 (Q), u 0 ∈ L 1 (R N ) are nonnegative, radially symmetric decreasing functions with respect to |x|. If for almost all t > 0, then the conclusion u # (|x|, t) ≺ v(|x|, t) still holds.
A priori estimates and best constants. We use the parabolic comparison results of [95] to obtain precise a priori estimates for the solutions of equation (55). One of these estimates is the so-called L 1 into L ∞ smoothing effect in the following form: The estimates were obtained in [51,52] (and in the non-power case in [94]). In the paper [96] we obtain the precise exponents and, what is important, the best constant C in the decay inequality. The calculation of best constants in functional inequalities is a topic of continuing interest in the theory of PDEs, both in the elliptic and evolution settings. Classical references to the calculation of best constants by symmetrization methods are Aubin and Talenti's computation of the best constants in the Sobolev inequality in [12,85]. Our calculation of a priori estimates with exact exponents and best constants is closely related to the sharp decay estimate for solutions of the porous medium/fast diffusion equation in [89,90]. When treating the linear case ϕ(u) = u, the estimates are called ultra-contractivity, see the book [49] where the importance of best constants is stressed for the applications in Physics. As a further application of the comparison techniques, optimal estimates with initial data in Marcinkiewicz spaces are obtained.
An important critical exponent appears repeatedly in the paper as a lower bound, (61) m c := (N − 2s)/N .
Since we are assuming m > 0 and 0 < s < 1, it does not appear in dimension N = 1 if s ≥ 1/2. Thus, we study the question of deciding the possible extinction of solutions in the range m < m c in terms of some norm of the initial data, and estimating the extinction time. First of all, we construct an explicit extinction solution of the fractional fast diffusion equation in this range of m's, formula (35). Then, we obtain optimal estimates by using comparison based on symmetrization. In this direction we improve significantly the results of the previous papers [27] and [52], by obtaining optimal estimates on the extinction time for data in Marcinkiewicz spaces. Finally, all the results are stable under the limit s → 1, where the standard diffusion case is recovered. See [52] for details on such limit.

KPP propagation and fractional diffusion
The problem goes back to the work of Kolmogorov, Petrovskii and Piskunov, see [68], that presents the most simple reaction-diffusion equation concerning the concentration u of a single substance in one spatial dimension, ∂ t u = Du xx + f (u). The choice f (u) = u(1 − u) yields Fisher's equation [57] that was originally used to describe the spreading of biological populations. The celebrated result says that the long-time behavior of any solution of u(x, t), with suitable data 0 ≤ u 0 (x) ≤ 1 that decay fast at infinity, resembles a traveling wave with a definite speed. In dimensions N ≥ 1, the problem becomes This case has been studied by Aronson and Weinberger in [9], where they prove the same result as the one-dimensional case. The result is formulated in terms of linear propagation of the level sets of the solution. In the case of the more general model the same result as before holds in the case of slow diffusion m > 1 (Vázquez and de Pablo [50]). Departing from these results, King and McCabe examined in [66] the case of fast diffusion m < 1 of equation (62). For (N − 2) + /N < m < 1, the authors showed that the problem does not admit traveling wave solutions by proving that level sets of the solutions of the initial-value problem with suitable initial data propagate exponentially fast in time.
On the other hand, and independently, Cabré and Roquejoffre ( [33]) studied the case of fractional linear diffusion u t (x, t) + (−∆) s u(x, t) = f (u), where (−∆) s is the Fractional Laplacian operator with s ∈ (0, 1) and concluded in the same vein that there is no traveling wave behavior as t → ∞, and indeed the level sets propagate exponentially fast in time. This came as a surprise since their problem deals with linear diffusion.
Motivated by these two examples of break of the asymptotic TW structure, we studied in [81] the case of a diffusion that is both fractional and nonlinear. More exactly, we consider the following reaction-diffusion problem (63) u t (x, t) + (−∆) s u m (x, t) = f (u) for x ∈ R N and t > 0, u(x, 0) = u 0 (x) for x ∈ R N , We are interested in studying the propagation properties of nonnegative and bounded solutions of this problem in the spirit of the Fisher-KPP theory. We assume that the reaction term f (u) satisfies where the exponent λ(N, s, m) is stated explicitly in the different ranges of the exponent m.
In our work we establish the negative result about traveling wave behaviour, more precisely, we prove that an exponential rate of propagation of level sets is true in all cases. Due to the nonlinearity, the solution of the diffusion problems involved in the proofs does not admit an integral representation as the case m = 1. Instead, we use as an essential tool the behavior of the fundamental solution of the Fractional Porous Medium Equation, also called Barenblatt solution, recently studied in [93]. This allows us to explain the propagation mechanism in simple terms: the exponential rate of propagation of the level sets of solutions (with initial data having a certain minimum decay for large |x|) is a consequence of the power-like decay behaviour of the fundamental solutions of the FPME. To be precise, the decay rate of the tail of these solutions as |x| → ∞ is the essential information we use to calculate the rates of expansion. This information is combined with more or less usual techniques of linearization and comparison with sub-and super-solutions. We also need accurate lower estimates for positive solutions of this latter equation, and a further self-similar analysis for the linear diffusion problem. The delicate details are explained in [81].
Main results. The existence of a unique mild solution of problem (63) follows by semigroup approach. The mild solution corresponding to an initial datum u 0 ∈ L 1 (R N ), 0 ≤ u 0 ≤ 1 is in fact a positive, bounded, strong solution with C 1,α regularity. Let us introduce some notations.
Here is the precise statement of our main results for the solutions of the generalized KPP problem (63).
Our main conclusion is that exponential propagation is shown to be the common occurrence, and the existence of traveling wave behavior is reduced to the classical KPP cases mentioned at the beginning of this discussion Equations similar to (63) present interest for various researcher groups, like the groups led by X. Cabré, J.M. Roquejoffre, H. Berestycki, F. Hamel, and P. Felmer.

Current work and comments
A number of related models, issues and perspectives on elliptic and parabolic equations involving fractional Laplacians and more general integral operators is under way.
• Numerics. The numerical analysis of the solutions of the FPME was started by Teso [86] where the case s = 1/2 was studied. The extension method that we use to implement the fractional Laplacian has a number of specific difficulties for s = 1/2 that we have addressed in a subsequent paper [87]. One of the main differences of our work is that we do not directly deal with the integral formulation of the fractional Laplacian; instead of this, we pass through the Caffarelli-Silvestre extension, mentioned above.
Previous works dealing with the numerical analysis of nonlocal equations of this type are due to Cifani, Jakobsen, and Karlsen in [46]. In particular, they formulate some convergent numerical methods for entropy and viscosity solutions. The numerical analysis of the elliptic PDE (−∆) s u = f in a bounded domain with zero boundary data via the extension method has been recently studied by Nochetto and collaborators using finite elements, [75].
• Disclaimer. We have not covered the extensive work on stationary states, i.e., elliptic equations of fractional type. Neither did we go into the probabilistic approach that comes from long ago and has been very actively pursued up to these days. Evolution equations involving the fractional Laplacian appear in many applied models and some of this directions are now actively pursued, like quasi-geostrophic flows, a topic that we have touched in the paper [53], but has not been addressed here. Another interesting direction concerns geometric flows with fractional operators. This list does not aim at being representative. Finally, we have not discussed the nonlocal diffusion evolution model treated in the monograph [5] which uses integral operators with kernels that do not resemble the fractional Laplacians and lead to a different theory.