Bifurcation diagrams and multiplicity for nonlocal elliptic equations modeling gravitating systems based on Fermi-Dirac statistics

This paper is devoted to multiplicity results of solutions to nonlocal elliptic equations modeling gravitating systems. By considering the case of Fermi-Dirac statistics as a singular perturbation of Maxwell-Boltzmann one, we are able to produce multiplicity results. Our method is based on cumulated mass densities and a logarithmic change of coordinates that allows us to describe the set of all solutions by a non-autonomous perturbation of an autonomous dynamical system. This has interesting consequences in terms of bifurcation diagrams, which are illustrated by a some numerical computations. More specifically, we study a model based on the Fermi function as well as a simplified one for which estimates are easier to establish. The main difficulty comes from the fact that the mass enters in the equation as a parameter which makes the whole problem non-local.

1. Setting of the problem. We consider a stationary solution of the drift-diffusion equation with a nonlinear diffusion based on some function H, coupled with the gravitational Poisson equation ∆φ = ρ on a bounded domain in R 3 , under appropriate boundary conditions. The simplest case, that we will call the (MB) case, corresponds to for Maxwell-Boltzmann statistics, yielding classical linear diffusion in (1). In this paper we shall study the nonlinear diffusion corresponding to Fermi-Dirac statistics, cf. [4,5,6,23]. In that case the function H is given by is a Fermi function. See Appendix A for more details. We shall refer to this case as the (FD) case. In order to deal with more explicit estimates, we shall also consider a simplified Fermi-Dirac model. This (sFD) case captures the asymptotic behavior of (FD) and is given by Connection between (FD) and (sFD) cases can be established when the positive parameters η and µ are such that We again refer to Appendix A for further details. The (FD) model was introduced by P.-H. Chavanis et al. in [11,12] in the context of astrophysical models of gaseous stars (also see [10] for a general review of the subject). It can be seen as a perturbation of the (MB) model, or linear diffusion model, and we shall prove that some features of the set of the stationary solutions are shared by the linear (MB) model and the nonlinear models (FD) and (sFD), if we take η > 0 sufficiently small (or µ > 0 sufficiently large).
We have several reasons to consider the (sFD) case: it has all qualitative features of the (FD) model, equivalents in the asymptotic regimes are much easier to control and numerically it avoids painful computations of Fermi functions, without loosing anything at the level of the mathematical results (qualitative behavior of the solutions) and their physical interpretation.
In order to use the cumulated mass technique, we shall assume that the domain under consideration is the unit ball B := {x ∈ R 3 : |x| < 1} and that the above equations are respectively supplemented with no-flux boundary conditions ∇H(ρ) + ∇φ · n = 0 on ∂B for the mass density (here n(x) = x/|x| for any x ∈ ∂B) and homogenous Dirichlet boundary conditions for the potential φ = 0 on ∂B .
We are interested in the stationary problem with a fixed mass constraint B ρ dx = M that can be solved by ρ = F (λ − φ) for some appropriate Lagrange multiplier λ. Since F in all considered cases is monotone increasing, the multiplier λ is uniquely defined. Here F = H −1 is the inverse of H, that is, in case (FD) of Fermi-Dirac statistics, and in case of (MB) Maxwell-Boltzmann statistics. In the (MB) case, the Lagrange multiplier is explicitly given by The function F has no simple expression in the (sFD) case. In all cases the stationary problem that we have to solve can be formulated as Let us start with a result in the (MB) case, which is easy to visualize on the bifurcation diagram expressing the dependence of the supremum norm of φ on the mass parameter. This result is a simple reformulation of a former result from [2]. The corresponding diagram exhibits a spiraling structure that can be seen in Fig. 1 (left) and accounts for the multiplicity of solutions which is reflected by the oscillating behavior of the branch in the bifurcation diagram: see Fig. 1 (right). Proposition 1. There exist two sequences (M n * ) n≥1 and (M * n ) n≥1 which are ordered and monotone: Sketch of the proof. We start by reducing (3), written with F (z) = e z to the autonomous, dynamical system and look for the branch of solutions s → (x(s), y(s)) such that lim s→−∞ (x(s), y(s)) = (0, 0), where x and y are defined, consistently with notation in (5) and (7) to be specified later, as x(log r) = 1 r r 0 s 2 e λ−φ(s) ds and y(log r) = r 2 e λ−φ(r) ∀ r > 0 .
Then there exists a unique heteroclinic orbit joining the points (0, 0) and (2, 2) as shown in Fig. 1. This heteroclinic orbit can be used to parametrize all solutions to (3). Since s → (x(s + s 0 ), y(s + s 0 )) is also a solution to (4), we may choose s 0 ∈ R such that 4π x(s 0 ) = M at least if M is in the admissible range 4π x(R) and thus get a solution with e λ = y(s 0 ). The parameters M n * and M * n correspond to the lower and upper values of M at which the bifurcation curve turns.
The proof of Proposition 1 is standard in the theory of gravitating systems. The reader is invited to check that all solutions are indeed described by the scheme and is invited to refer to [2] for more details. Numerically, the scheme allows one to compute all solutions by solving a simple ODE problem. It is also at the core of our approach for the nonlinear case and this is what we are now going to explain.
The main result of this paper is the following theorem. It is a multiplicity result for the Fermi-Dirac model (FD) and its simplified version (sFD) stating that for some values of the mass parameter we have at least the same number of solutions as for the Maxwell-Boltzmann case, if the parameter η > 0 is small enough. We consider the two sequences (M n * ) n≥1 and (M * n ) n≥1 that have been defined in Proposition 1 in the (MB) case.
if η > 0 is sufficiently small, there are at least 2n solutions of (3) in the (FD) and (sFD) cases.
In other words, the model corresponding to the Fermi-Dirac statistics or its simplified version at least partially inherits the spiraling structure of the model corresponding to the Maxwell-Boltzmann. See Figs. 1 (left) and 2. From the numerical computations it can easily be conjectured that a more precise description of the solution set can be achieved, with exact multiplicity. Hence it seems that for M ∈ (M n * , M n+1 * ) and for M ∈ (M * n+1 , M * n ) the exact multiplicity are respectively 4n + 1 and 4n − 1, in the asymptotic regime corresponding to η → 0 + . However, such a conjecture requires a by far more delicate analysis than the one we have done in this paper and is therefore still open. This can be summarized in the bifurcation diagrams for the simplified Fermi-Dirac and various values of η approaching 0. See Figs. 1 (right) and 3. More details on bifurcation diagrams and further numerical results will be given in Section 6.
For the convenience of the reader and to avoid further repetitions, let us synthetize our notation: (MB) Maxwell-Boltzmann statistics: where µ is given in terms of η by (2), f −1/2 has already been defined and Here the function R, which will be useful in the computations, is defined by R = 1/H , and we shall also consider P such that P (z) = z H (z). Unless specified log(1 + m) otherwise, η > 0 implicitly means that (FD) and (sFD) are under consideration, while statements corresponding to η = 0 apply to (MB). However, let us emphasize that (MB) is a singular limit as η → 0 + of the cases corresponding to η > 0. It is the main purpose of this paper to clarify this issue. For simplicity, we shall omit to mention the dependence on η and write R = R η only when emphasizing the dependence on η.
As we shall see, part of the branches depicted in the bifurcation diagrams converges as η → 0 + to the branches of (MB). For small values of the mass parameter both models share the same existence and uniqueness property of solution as was proved in [14] by the generalized Rellich-Pohožaev method. However, there is a major difference in the asymptotic behavior when comparing the Fermi-Dirac and the Maxwell-Boltzmann cases. In the (MB) framework, there is a critical mass above which no stationary solution exists. In the (FD) and (sFD) cases, no such critical mass appears as was shown in [18], and for any mass there exists a solution. Numerically, it is quite clear how the spiral of Fig. 1 gets regularized in Fig. 2.
It is also rather clear that our results can be extended at almost no cost to a large class of nonlinearities F depending on a parameter η, with appropriate properties, that converge to the exponential function in the limit as η → 0 + . From the physics point of view, however, it is the (FD) case that makes sense and this is why we have chosen not to cover the most general case but only the (FD) and (sFD) nonlinearities.
Before entering in the details, let us give a brief account of the literature on the subject. The reader is invited to refer to the references given in the quoted papers, especially for historical developments of the subject. In the three-dimensional (FD) case, it has been shown in [18] that there exists a global branch of solutions of (3) with arbitrary masses M > 0. This result is achieved by a variational approach as in [5,18,24] and by topological arguments (also see [22]) based on the fixed point is chosen in order to satisfy the mass constraint in (3). Still in the framework of Fermi-Dirac statistics, see [11,13,14,21,22] for stationary models and [7,8,19] for the corresponding equations of evolution. In the (MB) case, the problem reduces to the classical Gelfand problem, cf. [3,16]. A related family of problems can be considered when the ρ factor in (1) is replaced by ρ P (ρ), where the pressure function P is given in terms of H by P (ρ) = ρ H (ρ). Such equations are motivated by [11,12] and have been mathematically studied in [4,5,6,23]. This paper is organized as follows. In Section 2, we will generalize the change of variables that has been done in the proof of Proposition 1 for η = 0 to the case η ≥ 0 and get a non-autonomous dynamical system. In the next section we shall focus our attention on some a priori estimates for the dynamical system. Section 4 is devoted to the study of the dependence of the supremum norm of the density on the mass parameter. Finally, in Section 5 we establish some continuity results and prove the existence of multiple solutions (Theorem 2) as suggested by the bifurcation diagrams. Detailed numerical results and plots of the bifurcation diagrams in the (sFD) case are presented in Section 6, which also contains plots of quantities of physical interest like the free energy. Technical results concerning Fermi and related functions can be found in Appendix A.
2. The dynamical system. As a preliminary observation, we recall that, according to the symmetry result of Gidas, Ni and Nirenberg in [15,Theorem 1], any solution φ to (3) is radially symmetric when F corresponds to (FD), (sFD) or (MB). In this section we shall prove that all radial solutions to (3) can be parametrized by the solutions to some dynamical system. With a standard abuse of notation, we shall write a radial function of x as a function of |x|.
On the other hand, if we integrate the Poisson equation r −2 r 2 φ = ρ once, then by smoothness of φ we know that lim r→0 + r 2 φ (r) = 0 and thus get Altogether, we have obtained that and, using ρ = r −2 ζ (r) and R(z) := 1/H (z), we find that Following the computations of [2,20], we introduce the following change of variables and find by differentiation of x(s) = e −s ζ(e s ) that (6) can be rewritten as (e −2s y) e −s + e −s R e −2s y x = 0 , thus showing that the dynamical system s → (x(s), y(s)) obeys the equations y = 2 y − e 2s R e −2s y x .
If R(z) = z, which corresponds to (MB), we recover the autonomous dynamical system (4) that was obtained in the proof of Proposition 1. For (FD) and (sFD), the main difficulty is due to the fact that (8)-(9) is not autonomous. However, our strategy is built on the observation that R converges to the identity as η → 0 + .

3.
A priori estimates on the dynamical system. In this section we establish some a priori estimates on the solutions of (8)- (9). We start with some observations on invariant regions.
Lemma 4. Consider the dynamical system (8)- (9) with R corresponding to one of the three cases, (MB), (FD), or (sFD), for some η ≥ 0. Then R is continuous and such that 0 ≤ R(z) ≤ z for any z ∈ R + . As a consequence, the y = 0 axis is a stable manifold under the action of the flow, while the half-line y = 3 x, x > 0 is tangent to the unstable manifold. Moreover, all trajectories with y ≥ 3 x, x > 0 at s = 0 are out of the positive quadrant for some negative time or, to be more specific, are such that x(s) < 0 or y(s) < 0 for any s < 0 large enough.
Proof. The estimate R(z) ≤ z in the (FD) case will be shown in Lemma 9, in the Appendix. It is straightforward in the other cases.
that (3 x(t) − y(t)) e t < 0 for some t < 0, arbitrarily small, and again reach a contradiction by the estimate 3 x(s) e s ≤ 3 x(t) e t − y(t) (e t − e s ) , s < t if we assume that x(s) and y(s) are positive for any s < t. It should be noted that the crucial monotonicity of (3 x(t) − y(t)) e t for x(t) > 0 was used in the argument above to get the conclusion.
Notice that the existence of a branch of solutions follows from the existence theory for elliptic equations that can be found, among others, in [18] and [22]. 4. A priori estimates depending on the mass. Let us notice that readily follows from the definition of M = B ρ dx. Here B is the unit ball in R 3 .
Proof. Indeed, from (9) and H = 1/R we get the following equality Next we may integrate (9) and use lim s→−∞ y(s) = 0, according to Lemma 5, and the monotonicity of R (see the expression of the Fermi function f α in Appendix A in the (sFD) case), to get where the inequality follows from lim s→−∞ y(s) e −2s = ρ ∞ and y(s) e −2s ≤ 0, as shown in Section 3. Inserting this estimate into (15) we finally have which is exactly the claim of the lemma because of (10) together with R ≥ 0 implying monotonicity of H and Lemma 5 with its proof implying 3 m = 3 x(0) ≥ y(0).
Recall that solutions exist for arbitrary large masses in the (FD) and (sFD) cases, that is for η > 0, according to [18]. Estimates (13) and (14) provide an interesting qualitative property of the branch of solutions for large values of the mass, which is of interest by itself and provides some insight on numerical results of Section 6. The proof of the following corollary is straightforward by (13) and (14).
In order to emphasize the dependence on η ≥ 0 we shall denote the solution by (x η , y η ).
We know that lim t→−∞ A η (t) = lim t→−∞ B(t) = 0. From the above system we deduce that the estimates , hold for almost all t < 0. An integration on (−∞, t) shows that , for any t < 0. Altogether, using a Gronwall estimate, we finally arrive at which concludes the proof.
Proof of Theorem 2. Let ε > 0 and consider the enlarged singular set of masses 4π is smooth and locally invertible, if (x 0 , y 0 ) denotes the solution to (4) satisfying (17). This property is elementary (see the proof of Proposition 1 for details), and it is is also satisfied by the solution (x η , y η ) to (8)- Our numerical scheme goes as follows. We fix an ε > 0, small, and use the fact that y(s) is of the order of ε for s = t(ε) if e −2t(ε) ε ≈ p ∞ , which determines the dependence of t(ε) on ε. Hence our approximated solution is given by To emphasize the dependence on ε, we shall denote it by (p ε , q ε ). Figs. 4 shows the trajectories in the (p, q) phase space.  The bifurcation diagrams show that as the parameter η approaches zero, the solutions of simplified Fermi-Dirac (sFD) model parametrized by this number η approach the ones for the Maxwell-Boltzmann (MB) case, corresponding to the limit value η = 0, when the mass is in the admissible range for (MB). They also show that solutions with arbitrarily large masses exist as long as η is positive.
In practice, for numerical purposes, we have chosen ε = 10 −6 p ∞ . We may notice that our solutions come very close to (0, 0) for some s < t(ε) but are such that x(s) is negative for larger, negative values of s. The Gelfand problem is recovered for η = 0 (see for instance [14,16]), but Figure 5 clearly shows that a solution exists for any M > 0 as soon as η is positive. This corresponds well to the known existence results which either by variational or topological arguments yield the solution irrespectively of any value of the mass parameter. For the detailed analysis of this issue, see [18].
Solutions to (3) can be obtained by a fixed point method as in [18]. Alternatively, they can be characterized variationally. Hence, it should be noted that to get the minimal solution of (3), one can minimize the free energy functional of the solution (cf. [18]), namely  (1) is, at least formally, the gradient flow of F with respect to Wasserstein's distance according to the ideas introduced, e.g., in [17] and the reader is invited to check that F is monotone non-increasing along the flow defined by (1). This flows preserves the mass. Hence a minimizer of F under mass constraint is non-linearly dynamically stable.
Entropies for the isothermal model and for a model with fixed energy have been exhibited for instance in [6], in connection with many other papers dealing with generalized entropy (or free energy) functionals, see for instance [1,6,9]. Note that for the isothermal case, the entropies in [5] or [6] differ by the mass constant with the above one (which has no physical consequences since mass is conserved by the flow).
In the general case, we may notice that the entropy generating function β can be written as β(z) = z H(z) − P (z) , where the pressure P is such that P (z) = z H (z), while for the specific R corresponding to (sFD) case (see Appendix A for details) we have β(z) = z log z − z + 9 10 η z 5/3 .  In both cases, the bold curve corresponds to η = 0. which shares the same asymptotic as the function R in the (FD) case. The constants η and µ are related by (2) so that the function H(z) = log z + 3 2 η z 2/3 has the same behavior as in the case (FD) as z → ±∞.
Proof. A straightforward computation shows that we have z−R(z) = η z 5/3 1+η z 2/3 in the (sFD) case. The upper bound in the (FD) case follows by considering equivalents as z → ∞ with w = z η 3/2 and C = sup As for the lower bound, we observe that where the last line follows from an integration by parts. Coming back to the expression of R, R(z) ≤ z is equivalent to f −1/2 (t) ≤ 2 f 1/2 (t) with t = f −1 1/2 (2z/µ), which is precisely what we have just shown. c 2013 by the authors. This paper may be reproduced, in its entirety, for non-commercial purposes.