Some free boundary problems involving non-local diffusion and aggregation.

We report on recent progress in the study of evolution processes involving degenerate parabolic equations which may exhibit free boundaries. The equations we have selected follow two recent trends in diffusion theory: considering anomalous diffusion with long-range effects, which leads to fractional operators or other operators involving kernels with large tails; and the combination of diffusion and aggregation effects, leading to delicate long-term equilibria whose description is still incipient.

Introduction

The systematic mathematical study of free boundary problems is quite recent compared to the present importance of the field, and goes back to the last half of the past century. In this period, it has developed in many directions, which combine modelling of quite diverse natural phenomena, particularly in diffusion, wave mechanics and elasticity. A main topic in this broad landscape has always been the study of evolution processes where the state equations governing the different phases are nonlinear parabolic equations, possibly degenerate or singular at certain points or for certain values of the unknowns. This is the case of the well-known models like the Stefan problem for the common evolution of two immiscible fluids, the porous medium equation and other models to be mentioned later.

Free boundary models combine the difficulty of the analysis of the nonlinear PDEs with the difficulty in locating the separating interface between the phases, or free boundary in the absence of a theory of classical solutions that is generally unavailable.

So, in principle, we can solve the problem in some generalized sense but neither the state variables are guaranteed to be continuous nor the free boundaries to be reasonable hyper-surfaces. Even when the initial and boundary data are smooth, nonlinear processes involving free boundaries may develop singularities in finite time both in the solution and the free boundary. This is therefore a difficult problem for both analysis and geometry. In the last 50 years, considerable progress has been achieved and this is reported in different monographs, and in its recent trends in other articles of this issue. Several new models have been actively pursued in recent years with various levels of mathematical rigour versus practical applicability. In this article, we will discuss free boundary problems arising in two new scenarios, non-local diffusion and aggregation/swarming processes.

Non-local diffusion

The classical theory of diffusion is expressed mathematically by means of the heat equation, and more generally by parabolic equations, normally of the linear type; such approach 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. Mean curvature flows also belong to this broad category. Reaction diffusion has also attracted considerable attention.

In the last decade, there has been a surge of activity focused on the use of the 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 of 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 non-local character of the process is expressed.

More generally, research on non-local operators has a tradition related to stochastic processes and has witnessed a rapid expansion in the last decade when it has attracted the attention of PDE experts who brought new problems and techniques into the field. The area poses new challenges to both pure and applied mathematicians and is now at the interface of at least three wide fields, functional analysis, stochastic processes and partial differential equations, while providing a new paradigm in scientific modelling. This interaction between mathematics and applications is giving rise to new concepts and methods and is expected to produce new challenging mathematical problems for many years to come.

Fractional operators

Though there is a wide class of interesting non-local operators under scrutiny, both in theory and applications, a substantial part of the current work deals with diffusion modelled 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 |ξ|2. Thus, if g is a function in the Schwartz class in , N≥1, we write (−Δ)g=h if 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 0C is a normalization constant, with precise value that is found in the literature. In the limits s→0 and s→1, it is possible to recover, respectively, the identity or the standard minus Laplacian, −Δ (cf. [1]). Remarkably, the latter one cannot be represented by a non-local formula of the type (2.1). It is also useful to recall that the operators (−Δ)−, 02] and Stein [3]. 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. These and other definitions are equivalent when dealing with the Laplacian on the whole space . The correct definitions for the operators defined on bounded domains admit several options that are being investigated at the moment [4,5]. The interest in these fractional operators has a long history in probability and other applied sciences. 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; see for instance [6-10], where further references can be found.

A part of the current research concerns linear or quasi-linear equations of elliptic type. This is a huge subject with well-known classical references. Even if it includes free boundary problems, mainly of the obstacle type [11] it will not be discussed here in itself for lack of space, because we want to present free boundaries in evolution problems.

Linear evolution processes of anomalous type

We will study evolution models that arise as variants of the heat equation paradigm, including however long-range effects. Thus, the difference between the standard and the fractional Laplacian consists in taking into account long-range interactions (jumps) instead of the usual interaction driven by close neighbours. The change of model explains characteristic new features of great importance, like enhanced propagation with the appearance of fat tails at long distances: such tails greatly differ from the typical exponentially small tails of the standard diffusion, and even more from the compactly supported solutions of porous medium flows. Moreover, the space scale of the propagation of the distribution is not proportional to t1/2 as in Brownian motion, but to another power of time, that can be adjusted in the model; this is known as anomalous diffusion. Anomalous diffusion is nowadays intensively studied, both theoretically and experimentally, since it 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 (e.g. [12-14]). The fractional Laplacian operators of the form (−Δ), σ∈(0,2), are actually the infinitesimal generators of stable Lévy processes [12,15,16]. The standard linear evolution equation involving fractional diffusion is which is the main model for anomalous diffusion. The equation is solved with the aid of well-known tools, like Fourier transform. Posed in the whole space, it generates a semigroup of ordered contractions in and has the integral representation where K has Fourier transform . This means that, for 017]. When , F is explicit If s=1 the function K is the Gaussian heat kernel, which has a negative square exponential tail, i.e. a completely different long-distance behaviour.

From the point of view of the present investigation, this fractional linear model and other variants, though useful and heavily studied, fail to meet the requirement that we are looking for, that is, generating free boundaries in evolution. On the contrary, the solutions have large densities at long distances.

Nonlinear and non-local 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 have been proposed as nonlinear counterparts of the linear fractional heat equation, and combine Laplace operators and nonlinearities in different ways. Let us mention the two popular models in the recent PDE literature, and how they relate to our stated goal.

— Model I. A quite natural option is to consider the equation with 00. This is mathematically the simplest fractional version of the standard porous medium equation (PME) ∂u=Δ(u), that is recovered as the limit s→1 and has been extensively studied (cf. [18,19]). We will call equation (3.1) the fractional porous medium equation (FPME) as proposed by de Pablo et al. [20,21]. From the point of view of our present study, this model is quite interesting because it offers a possible balance between a porous medium nonlinearity, which for m>1 implies finite propagation and free boundaries, and a fractional-type diffusion where the speed of propagation is infinite and fat tails arise. The question is: which effect wins? In [20,21] we have answered the question in the sense of infinite propagation and no free boundaries (we are considering non-negative solutions of the evolution problem).

Nowadays, a reasonably complete theory of nonlinear diffusion with long-range effects has been developed for this model, see for instance [22-25], or similar models like [26], but we have to look elsewhere for free boundaries.

— Model II. A second model of nonlinear diffusion with fractional operators has been studied by the second author in collaboration with Caffarelli and does give rise to the occurrence of finite propagation and free boundaries. This alternative model is derived in a more classical way from the PME since it is based on the usual Darcy law (i.e. the velocity of the particles is assumed to be the gradient of a pressure function, v=−∇p) with the novelty that the pressure is related to the density by an inverse fractional Laplacian operator, . Putting these two facts into the continuity equation, ∂u+∇⋅(uv)=0, the model takes the form a nonlinear fractional diffusion equation of porous medium type where is the Riesz operator that typically expresses the inverse to the fractional Laplacian, . This has been studied by Caffarelli & Vázquez [27,28]. In dimension N=1, this model was studied by Biler et al. [29,30] as a model for the propagation of dislocations as proposed by Head [31]. Since it was extensively reported in a survey paper [9] we will only point out the established theory of weak solutions, and the lack of uniqueness results in several dimensions, and concentrate on the existence and properties of the free boundaries. Self-similar solutions exist and their existence and properties are very illustrative of the phenomenon of finite propagation.

Finite propagation for Model II. Solutions with compact support

One of the most important features of the PME and other related degenerate parabolic equations is the property of finite propagation, whereby compactly supported initial data u0(x) give rise to solutions u(x,t) that have the same property for all positive times, i.e. the support of u(⋅,t) is contained in a ball B(0) for all t>0. One possible proof in the case of the PME is by constructing explicit weak solutions exhibiting that property (i.e. having a free boundary) and then using the comparison principle, that holds for that equation. Since we do not have such a general principle here, we have to devise a valid form of comparison; we show that comparison works if we use a special family of excessive supersolutions technically called ‘true supersolutions’. The technique has to take into account the peculiar form of the integral kernels.

We begin with N=1 for simplicity. We assume that our solution u(x,t)≥0 has bounded initial data u0(x)=u(x,t0)≤M with compact support and is such that u0 is below the parabola a(x−b)2, a,b>0, with graphs strictly separated. We may assume that u0 is located under the left branch of the parabola. We take as comparison function which is a travelling wave moving to the right with speed C that will be taken big enough. Then we argue at the first point and time where u(x,t) touches the left branch of the parabola U from below. The key point is that if C is large enough such contact cannot exist. The formal idea is to write the equation as and observe that at the contact we have u≥U=2aC(Ct−x+b), while u=U=−2a(Ct−x+b), so the first can be made much bigger than the second by increasing C. The influence of p and p as well as u is controlled, and then we conclude that the equation cannot hold if C is large enough. The argument can be translated for several dimensions. Here are the detailed results proved in [27].

Theorem 3.1

Let and assume that u is a bounded solution of equation (3.2) with 0≤u(x,t)≤L, and u0 lies below a function of the form U0(x)= A e−, A,a>0. If A is large then there is a constant C>0 that depends only on (N,s,a,L,A) such that for any T>0 we will have the comparison If a similar statement is true but C is not a constant but some increasing function of time.

A similar finite propagation result is true in several space dimensions. The study of the corresponding free boundaries is an open topic. Let us mention in passing that much effort has been devoted in [27,32,33] to prove the basic sequence of regularity results, according to which initial data in L1 give rise to solutions that belong to L, for all positive time. Moreover, there exists a positive constant C such that for every t>0 with α=N/(N+2−2s), γ=(2−2s)/(N+2−2s). Finally, it is proved that bounded weak solutions u≥0 of the initial value problem are uniformly continuous on bounded sets of s<1. Indeed, they are C continuous with a uniform modulus.

Self-similar solutions for Model II

Next we want to study the large time behaviour of this model following [28], since it sheds some light on the way finite propagation works. The first step is constructing the self-similar solutions that will serve as attractors. Surprising result: their density will be compactly supported, while their pressure will be positive everywhere with the typical fat tails far away.

Rescaling for the fractional porous medium equation

Inspired by the asymptotics of the standard porous medium equation, we define the rescaled (also called renormalized) flow through the transformation with new time . We also put y=x/(t+1) as rescaled space variable. In order to cancel the factors including t explicitly, we get the condition on the exponents α+(2−2s)β=1. Here we use the homogeneity of in the form . From physical considerations, we also impose the law that states conservation of (finite) mass, which amounts to the condition α=Nβ, and in this way we arrive at the precise value for the exponents: We also arrive at the nonlinear, non-local Fokker–Planck equation which is the equation for the renormalized flow. In all the above calculations, the factor (t+1) can be replaced by t+t0 for any t0>0, or even by plain t.

Stationary renormalized solutions

It is important to concentrate on the stationary states of the new equation, i.e. on the solutions V (y) of where a=β/2, and β is defined just above. Since we are looking for asymptotic profiles of the standard solutions of the FPME we also want V ≥0 and integrable. The simplest possibility is integrating once to get The first equation gives an alternative choice that reminds us of the complementary formulation of the obstacle problems. Indeed, if we solve the obstacle problem with fractional Laplacian, we will obtain a unique solution P(y) of the problem: with 0C−a|y|2, where C is any positive constant and a=β/2. Note that −Δf=2Na=α. For uniqueness, we also need the condition P→0 as . Fortunately, the corresponding theory has been developed by Caffarelli and collaborators (cf. [34,35]). The solution is unique and belongs to the space H− with pressure in H. Moreover, it is shown that the solutions have P∈C1, and V ∈C1−.

Note that for C≤0 the solution is trivial, P=0, V =0, hence we choose C>0. We also note the pressure is defined but for a constant, so that we could maybe take as pressure instead of P so that ; but this does not simplify things since P→0 implies that as . Keeping thus the original proposal, we get a one parameter family of stationary profiles that we denote V (y). These solutions of the obstacle problem produce correct weak solutions of the FPME with initial data a multiple of the Dirac delta for the density, in the form It is what we can call the source-type or Barenblatt solution for this problem, which is a profile V ≥0. It is positive in the contact set of the obstacle problem, which has the form , and is zero outside, hence it has compact support.

On the other hand, the rescaled pressure P(|y|) is always positive and decays to zero as according to fractional potential theory [3]. The rate of decay of P as turns out to be P=O(|y|2), a fat tail.

Exact calculation of density profiles

Biler et al. [29] studied the existence and stability of self-similar solutions in one space dimension. Recently, Biler et al. [30] obtained the explicit formula for a multi-dimensional self-similar solution in the form with α=N/(N+2−2s) as before. The derivation uses an important identity for fractional Laplacians which is found by Getoor [36]: (−Δ)(1−y2)+=K if σ∈(0,2]. Here we must take σ=2(1−s). According to our previous calculations ΔP=−α on the coincidence set, hence c1=α/K. Let us work a bit more: using a similar scaling to (3.3), we arrive at the following one-parameter family of self-similar solutions: where and C1>0 is a free parameter that can be fixed in terms of the mass of the solution . This is the family of densities that corresponds to the pressures obtained above as solutions of the obstacle problem. Let us finally mention that explicit formulae for self-similar solutions to related fractional equations have also recently been obtained [37].

Asymptotic behaviour for Model II

The next step is to prove that these profiles are attractors for the rescaled flow. The study is done in terms of the rescaled flow that is much more workable than the original evolution. The following result is proved in [28].

Theorem 3.2

Let u(x,t)≥0 be a weak solution of the Cauchy problem for equation (3.2) with bounded and integrable initial data such that u0≥0 has finite entropy, i.e. Let v(y,τ) be the corresponding rescaled solution to (3.4). As we have The constant C is determined by the rule of mass equality: . In terms of function u, this translates into the convergence as

The constants α and β are the self-similar exponents defined before. The proof uses entropy dissipation methods that adapt perfectly to the problem. Sharper results on this convergence are given in [38] where more references can be found. See also below in this respect.

Extension of the nonlinear fractional models

(1) There is way of generalizing the two previous Models I and II, to accept two exponents m and p: so that the comparison of both models happens on symmetric terms. The question we want to solve is deciding between finite and infinite propagation in terms of the exponents m and p. Progress on this issue is reported in [39,40].

(2) Limits. The limit of Model II when s→1 is quite interesting since one obtains a variant of the equation for the evolution of vortices in superconductivity studied in [41-44]. The understanding of this limit has been done by the second author in collaboration with Serfaty [45] and is related to work by Bertozzi et al. on aggregation models [46,47]. This is a good point to connect to the next issue. On the other hand, the limit of Model I leads to a free boundary problem of the mesa type reported in [23,48].

Aggregation and swarming

The collective behaviour of individuals in animal grouping, also called swarming, has recently been modelled by mean field PDEs (see [49] and references therein). They consist of continuum descriptions of individual based models (IBMs; see [50-52] for applications to fish schools and software simulation of swarms). The motion of individuals in IBMs is governed by systems of ODEs imposing an asymptotic speed, cruise speed, for particles. One of the most celebrated models of this kind was introduced in [53]. The authors propose to introduce a social force between individuals based on the empirical observation that there are two basic mechanisms of interaction: an inner repulsion zone to keep a comfort area or to avoid collisions and an outer attraction area since individuals do want to socialize and be close to each other. These two tendencies are modelled via a pairwise potential directed along the position between the particles and depending on the interparticle distance. Therefore, this model consists of where are the positions and velocities of the individual particles and α,β are effective values for self-propulsion and friction forces (see [53] for more discussions). These models lead to kinetic equations in the limit of large number of particles as derived in [54]. Some of the patterns that the simulation of these IBMs produce are the so-called flock patterns. The flock solution corresponds to the case in which all particles move with the same velocity vector v=v0 with |v0|2=α/β while the positions satisfy that the forces are balanced, that is, The flocking solution corresponds then to the translational movement of the shape formed by the positions with velocity vector v0. When the number of particles tends to infinity, these discrete flocking solutions are approximating continuum density profiles satisfying at least over the support of . In fact, the connection between the second-order microscopic model (4.1) and the first-order model given by has been clarified in [55,56]. In these articles, it is shown that the stability of flock solutions for (4.1) is equivalent to the stability of the spatial shape of the flock as steady solution to (4.2). A similar result at the level of the continuum equations is lacking. However, it makes sense to find all the spatial configurations leading to flocks or equivalently all stationary solutions to (4.2). Among those, the most stable configurations are of special interest, in view of the connection through the stability of the solutions at least at the discrete level. The system (4.2) has a natural energy that is dissipated along its evolution given by The continuum version of this energy, suitably scaled, leads to where μ is a probability measure. This energy is dissipated along the evolution of the continuum mean-field limit of the system of ODEs in (4.2), called the aggregation equation where ρ is the density function of the measure μ. Therefore, minimizers of the energy (4.3) on the set of probability measures should be among the most stable configurations. Finally, let us comment that both the non-local diffusion equation (3.2) and its rescaled version (3.4) are particular instances of the aggregation equation (4.4) for suitably chosen potentials W. We will see in the rest of this section a summary of the results concerning the minimizers of the interaction energy (4.3), especially in its relation to free boundary problems, particularly to the obstacle problem.

Connection to optimal transport

We start by noting that the discrete system (4.2) has the structure of a gradient flow with respect to the Euclidean distance. In fact, this structure is shared by the continuum equation (4.4) in a suitable sense. We need to be precise about the definition of distance between probability measures. We recall that is the set of Borel probability measures on and we denote by the family of Borel subsets of . The support of a measure is the closed set defined by A family of distances between probability measures has been classically introduced by means of optimal transport theory. We will review briefly some of these concepts; we refer to [57] for further details. A probability measure π on the product space is said to be a transference plan between and if for all . If , then denotes the set of admissible transference plans between μ and ν. Informally, if π∈Π(μ,ν) then dπ(x,y) measures the amount of mass transferred from location x to location y. We recall that the Euclidean Wasserstein distance d2 between two measures μ and ν is defined by Note that for the set of probability measures with finite moments of order 2. It is classical that equation (4.4) is a gradient flow of the energy (4.3) in the space of measures endowed with the distance d2 [58,59]. This property is essential to obtain many qualitative properties of the solutions and to deal with possible singular interaction potentials (see [46,47,60-66] and references therein).

In our context of minimizing the energy functional (4.3), this means that we look for local minimizers in d2. (Local) Minimizers of (4.3) should correspond to equilibrium configurations for the evolution equation obtained by steepest descent of the energy. However, being a functional on probability measures, the steepest descent has to be understood in the Wasserstein sense as in [67] by writing (4.4) as This evolution equation can make sense in the set of probability measures depending on the regularity of the potential W. However, if the potential is singular at the origin the well-posedness only happens in L spaces (see [60,68,69] for more details). For steady configurations, we expect ∇ψ=0 on the support of μ, owing to the formal energy dissipation identity for solutions, i.e. where μ is any solution at time t of (4.6) and ψ=W*μ its associated potential. Therefore, the points on the support of a local minimizer μ of the energy E should correspond to critical points of its associated potential ψ. This fact is made rigorous in the following proposition.

Proposition 4.1 ([64], theorem 4)

Assume that W is a non-negative lower semi-continuous function in . If μ is a d2-local minimizer of the energy, then the potential ψ satisfies the Euler–Lagrange conditions given by

(i) ψ(x)=(W*μ)(x)=2E[μ] μ-a.e.

(ii) ψ(x)=(W*μ)(x)≤2E[μ] for all x∈supp(μ).

(iii) ψ(x)=(W*μ)(x)≥2E[μ] for a.e. .

These conditions simplify to if μ is absolutely continuous with respect to the Lebesgue measure. These Euler–Lagrange conditions can also be interpreted in terms of game theory [70-72]. In fact, they imply that the potential ψ achieves its global minima on the support of μ, and therefore μ is also characterized as In game theory, the probability measure μ is interpreted as the strategy chosen by the agent x, the potential is interpreted as encoding the interactions between the agents, and the minimization is thought about the optimal condition (the game) that the agents are looking for. Therefore, the Euler–Lagrange conditions are interpreted in game theory as the lack of interest of any player in the game to a change of strategy leading to the concept of Nash equilibria. As a summary, local minimizers of the potential energy are among Nash equilibria within the game theoretical viewpoint.

Let us concentrate on the quest of local minimizers for repulsive/attractive potentials of the form with a>b>−N to be repulsive at the origin and attractive at infinity potential which is locally integrable. We allow a or b to be zero with the understanding that .

As mentioned before, the case a=2 and b=2−N, black point in figure 1, is of particular interest. It corresponds to Newtonian repulsion and it has been repetitively rediscovered that the unique, up to translations, global minimizer of the interaction energy is the characteristic function of a Euclidean ball with suitable radius [73-75]. The uniqueness, up to translations, of the global minimizer for more singular than Newtonian repulsion, a=2, b=2s−N with 028] via the connection to a classical obstacle problem; see also [75]. This strategy was also used in [45] to treat again the case s=1 for the evolution problem as in [47]. If the parameters are in the range of the green line in figure 1, i.e. a>2 or 2−N76,77]. Summing up, the solution support is limited by a free boundary, whose location (a sphere) is to be determined as part of the problem.

Figure 1.

Parameters for minimizers of the interaction energy with potential of the form (4.8) given by densities. See the text for a full explanation.

Moreover, the main result in [78] shows that in the range a>0 with b=2s−N with 0

Let us conclude this part by mentioning that whenever we can fully characterize the unique, up to translations, global minimizers of (4.3), they are the long-time asymptotics for the aggregation equation (4.4); see [27,28,45,47] for the known different cases (also §3). These cases are essentially reduced to a=2 and b=2s−N with 038,79]. Finally, let us remark that finding sufficient conditions on the interaction potential W to ensure the existence of compactly supported global minimizers is a very interesting question that has been recently solved in [80]; see also [81-83]. However, we do not know yet about their properties in terms of asymptotic stability for the aggregation equation (4.4).

More on the connection to free boundary problems and the obstacle problem

We will now focus on a topic already introduced in §3b. A classical obstacle problem [34] for a differential operator consists of finding a function such that where are just given functions. Here f is the obstacle to Φ, and thus whenever Φ does not touch the obstacle, then the first condition implies that Φ should be a solution to the PDE: . To make the connection to the Euler–Lagrange condition clearer, let us take the particular case of the classical potential encountered in semiconductors [84,85] given by with N≥1. Here, c is the constant such that the first part of the potential is the fundamental solution of the Laplacian operator Δ in N≥1. Assume that is a d2-local minimizer of the energy associated with this potential, then the potential associated to ρ, given by ψ=W*ρ is well defined as a function that satisfies in the sense of distributions and thus Δψ≥N in and Δψ=N outside the support of ρ. In conclusion, the potential ψ associated with the local minimizer satisfies a classical obstacle problem in with operator , g=N, and f=2E[ρ] due to (4.7). Let us also remark that the obstacle problem (3.6) associated with the scaled version of the fractional diffusion equation (3.4) can also be recast in the form (4.9). We will do this next but in more generality.

The basic idea for the Laplacian can be generalized and used to get properties on the potentials associated with local minimizers of a larger family of interaction potentials W without resorting to the local minimizers themselves. Let us review quickly the main strategy in [78]. The key fact is that if the repulsive strength at zero of the potential is given by |x|/b, b=2s−N with 078]. This allows us to deduce that for any point x0∈supp(μ) holds. Here, ε is the size of the ball in d2 distance where μ is a local minimizer. Moreover, for the points in the support, we also get Since |x|/b in this range is the fundamental solution, up to constants, of the fractional diffusion operator , we can write in the distributional sense that In particular, since μ is a non-negative measure, we deduce Furthermore, if x∈B(x0) is such that ψ(x)>ψ(x0), (4.11) implies that x∉supp(μ). We deduce Collecting (4.10) and (4.11), we have, at least formally, the following result.

Theorem 4.2 ([78])

For all x0∈supp(μ), the potential function ψ associated with a d2-local minimizer of the energy (4.3) is equal, in B(x0), to the unique solution of the obstacle problem where C0=ψ(x0) and F(x)=−(−Δ)(|x|/a)*μ. Furthermore, we can recover the local minimizer of the energy μ as μ=−Δψ+F.

Most of the results mentioned above follow this strategy to make use of the regularity and existence theorems for obstacle problems available in the literature [34,35,86-88]. Let us remark that in the case of a power-law potential (4.8) with a=2 and 028]. As before, theorem 4.2 is true for more general potentials behaving as the repulsive power |x| at the origin with suitable growth conditions at infinity [78].

It is an open question how to use obstacle problem techniques in cases in which the repulsion is not too singular, for instance for power-law potentials when (4.8) with a=2 and b=2s−N with now s>1. The main problem we face is that the inverse operator of the convolution is no longer as nice as it was for fractional diffusion. However, we are currently investigating the case when s=1+ϵ with 0<ϵ<1, i.e. in the range corresponding to the blue line in figure 1. In this particular case, we consider an intermediate obstacle problem by defining Ψ to be the solution to which again has nice regularity properties since 0<ϵ<1. Now, we define our potential ψ as the Newtonian potential associated with Ψ, that is, we solve Ψ=−Δψ. The idea is that by integrating twice the solution of the obstacle problem (4.13), we get a candidate to be a solution to (4.12).