On the Effects of Bohm's Potential on a Stationary Macroscopic System of Self-Interacting Particles

We consider a macroscopic model describing a system of self-gravitating particles. We study the existence and uniqueness of non-negative stationary solutions and allude the differences to results obtained from classical gravitational models. The problem is considered on a bounded domain up to three space dimension, subject to Neumann boundary condition for the particle density, and Dirichlet boundary condition for the self-interacting potential. Finally, we show numerical simulations that affirm our findings.

1. Introduction. Consider a stationary macroscopic system of self-interacting particles with Bohm's porential, which describe the normalized density n ≥ 0, the quasi Fermi-level F , −div (n∇F ) = 0 in Ω, n ∂ ν F = 0 on Γ, (1b) and the potential Φ due to self-interaction in a particle system, where ν is the outer normal to the convex, bounded domain Ω ⊂ R d , d ≤ 3 with Lipschitz boundary Γ, > 0 is the scaled Planck constant and |σ| ∈ [0, ∞) is the mass of the system of self-interacting particles, where sign(σ) dictates the nature of the interaction involved. In this case, positive mass σ > 0 would indicate the presence of self-attraction, while negative mass σ < 0 indicates self-repulsion. By passing to the limit → 0, we formally recover either the classical drift-diffusion equations (σ < 0) [1] or a model for a system of self-gravitating particles (σ > 0) [3].
A simple observation of (1b) suggests that for any positive density n > 0, we obtain a solution F ∈ R. In fact, any constant function is a solution. However, we shall see below that this constant solution is fixed due to the normality of n.
In order to bring this system of equations into a system similar to that of the classical equations for self-interacting particles (c.f. [1,3]), we introduce the quasi potential u := log n − F . Assuming n > 0, we insert F = log n − u into (1b) to obtain This allows us to relate n and u via n = αe u for some constant α > 0. Since n is normalized, i.e., Ω n dx = 1, we deduce that α = e u −1 L 1 (Ω) . Notice that, by fixing α, we also fix the quasi Fermi-level, which is explicitly given by F = log(α).
Introducing this into (1) leads to the coupled system for (u, Φ) given by an elliptic equation with natural gradient growth for u ∆u + u = 2 4 |∇u| 2 + σΦ in Ω, ∂ ν u = 0 on Γ, (3a) and the equation for the potential Φ, Note that system (3) is equivalent to system (1) if n > 0, or equivalently, if u is an essentially bounded function. The existence of bounded weak solutions (u, Φ) will be shown for (3), thereby implying the existence of solutions for (1).
We introduce the short hand X to denote the space Theorem 1. Let d ∈ {2, 3} and the mass |σ| ∈ [0, +∞) be given. Then problem is a solution of (1). Moreover, there exists θ ∈ (0, 1), such that θ ≤ n ≤ 1/θ. Observe that σ > 0 can be chosen arbitrarily large as opposed to classical selfgravitating particles, where a threshold for existence exists. In this sense, system (1) can be thought of as a regularization of the classical self-gravitating system.
Let us discuss the techniques used to show existence of solutions for system (3). The solvability of (3a) and its variants with homogeneous Dirichlet boundary data were shown in the papers [7,8,12]. Moreover, if Φ ∈ L p (Ω) with p > d/2, then u ∈ X . Adopting the methods used in [7], we show in Section 2 that this holds true also for homogeneous Neumann boundary data. In this case, systems (1) and (3) are equivalent. Furthermore, we obtain from [6] the following sharp estimates for functions in the space W 2,1 ∆,0 = D(Ω) ∆· L 1 (Ω) .
where ω d−1 is the measure of the unit sphere in R d . The constants given above are the best possible, independently of the domain.
Here we used L exp (Ω) to denote the Zygmund space, whose elements f satisfy Ω e λf dx < ∞ for some λ = λ(f ) > 0, and L p,∞ (Ω) to denote the classical weak-L p space (see [2]). These spaces are related to each other and to the L p spaces by the following continuous embeddings for 1 < p < ∞, The estimates in Proposition 2 necessarily implies that the L q -norm of Φ, for any q ∈ (1, d d−2 ), is uniformly bounded for d ≥ 2. In particular, we have where c q is the embedding constant obtained from (4).
Therefore, an application of the Schauder fixed point theorem on a self-mapping for the potential Φ leads to its existence in L q (Ω) for some q ∈ ( d 2 , d d−2 ), and consequently also for u and n, which yields the existence result. For sufficiently σ we further obtain uniqueness of solutions, given by the following result.
To our knowledge, this estimate for uniqueness appears to be new for both attractive and repulsive potentials. For σ < 0, it is known that in the classical case = 0, uniqueness depends on the smallness of the applied voltage. As a matter of fact, the performance of many semi-conductor devices (thyristors) depends on the existence of multiple solutions (c.f. [14] and references therein).
Unfortunately, this estimate is not sharp, since for d = 2 and σ > 0, it is known that uniqueness is valid for σ < µ 2 (c.f. [5]). Nevertheless, the above estimate provides a convenient relationship between uniqueness and the Bohm potential ( > 0).

Elliptic equation with natural gradient growth.
We begin by providing results for the subproblem (3a) required to prove Theorem 1.
Definition 4. A function u ∈ X is said to be a solution of (3a) if it satisfies for every function ϕ ∈ X .
Before proving the theorem, we state two results regarding the regularity of the solution u, whose proofs can be found in the appendix.
Proof of Theorem 5. Define φ: R → R by φ(s) = (e |s| − 1) sign(s) and the cut-off function T k : R → R by T k (s) = max{−k, min{s, k}} for some k ∈ R. We begin by considering the auxiliary problem for u k ∈ X : Since the right-hand side is bounded, the existence of a bounded solution for (9), k ∈ N, may be deduced from classical results (see for example [11] for the existence and [10] for the boundedness). Due to the fact that T k (|∇u k | 2 ) ≤ |∇u k | 2 and |T k (Φ)| ≤ |Φ| along with Lemma 6, there exists a function u ∈ X such that In order to pass to the limit in (9), we still need to show that ∇u k → ∇u in L 2 (Ω), i.e., the strong convergence of the gradients of u k in L 2 (Ω). To do so, we test (9) with For the first term on the left-hand side we have As for the second term on the left hand-side, we have For the first term on the right-hand side we have Due to the compact embedding H 1 (Ω) → L 2 (Ω), we get a convergent subsequence, denoted again by {u k }, such that u k → u in L 2 (Ω). Consequently we obtain yet another subsequence, denoted again by {u k }, such that u k (x) → u(x) for a.e. x ∈ Ω, which implies the almost everywhere convergences From Lebesgue's dominated convergence for these sequences and their boundedness in L ∞ (Ω), we have the strong convergences Passing to the limit in (10) yields which necessarily implies that ∇u n → ∇u in L 2 (Ω). Therefore, passing to the limit in (9) yields the solution u ∈ X satisfying (6). The fact that e λu ∈ X for every λ ≥ 0 follows directly from Lemma 6.
3. Proof of Theorem 1. As drafted out above, we use the Schauder fixed point theorem (c.f. [9,Corollary 11.2]) to facilitate the proof. We define the closed, convex and bounded subset of L q (Ω) and with µ d as given in (5). For a given w ∈ M , we consider the auxiliary problem for Φ given by This induces a compact mapping simply due to the continuity of the solution operators between their respective spaces and the compact embedding H 1 (Ω) → L q (Ω). Indeed, since w ∈ L q (Ω), we obtain a solution u ∈ X of (11a) as a result of Theorem 5. Inserting u into (11b) and solving for Φ gives us Φ ∈ X with Φ ≥ 0 due to standard theory of elliptic equations. Furthermore, we have H: M → M simply due to the estimates in (5). A direct application of the Schauder fixed point theorem concludes the proof.
4. Proof of Theorem 3. Let A: X → H 1 (Ω) * denote the operator defined by Then the weak formulation of (1a) can be written as n ∈ X : With similar arguments to those by Pinnau, Unterreiter in [14, Theorem 26] the operator A is well-defined. Moreover, for a fixed ϕ ∈ H 1 (Ω), the Gateâux derivative of A(·), ϕ : X → R at a point n ∈ X in any direction h ∈ X exists and is given by , be two solutions of (1) and set the difference to be (δn, δF, δΦ) : Setting ϕ = δn we obtain by subtraction from (12) A(n 1 ) − A(n 2 ), δn = σ δΦ, δn + δF, δn = σ δΦ, δn , where we used the fact that δF ∈ R and Ω δn dx = 0.
Let n τ = n 1 − τ δn, τ ∈ [0, 1] be the convex combination of n 1 and n 2 . Since the function [0, 1] τ → A(n τ ), δn ∈ R is differentiable, we obtain, by the mean value theorem, a τ ∈ (0, 1) such that Now, in order to estimate the first term on the left-hand side from below, we use a variant of the result obtain in [14,Lemma 24].
For the second term on the left-hand side, we apply Hölder's inequality to obtain where we used the fact that Ω n τ dx = 1. Altogether we have the estimate for some s ∈ [1, ∞) satisfying the requirements of Proposition 7. Note that by subtraction, δΦ ∈ X solves the problem −∆δΦ = δn in Ω, δΦ = 0 on Γ.
As before, the estimates in Proposition 2 yield for q ∈ ( 4d d+2 , d d−2 ) the inequality Due to the constraint placed on q, we have the Sobolev embedding H 1 (Ω) → L q (Ω), where q = q q−1 , thereby allowing us to choose s = q in (14). Putting together all the inequalities obtained above, we have Applying Young's inequality we obtain Using the continuous embedding L q (Ω) → L 1 (Ω), we finally arrive at In conclusion, for sufficiently small mass we obtain uniqueness for n, and consequently for Φ.

The semi-classical limit.
This result is well-known for the case σ < 0, i.e., for the quantum drift-diffusion equations (c.f. [1]). Therefore, we restrict ourselves to the case σ > 0. Furthermore, we consider only the case d = 2.
According to [4] the free energy functional where Φ is a solution of the Poisson problem (1c) attains a minimum n 0 in the set P = n ∈ L 1 (Ω) | n ≥ 0, Ω n dx = 1, n log n ∈ L 1 (Ω) .
Moreover, the minimizer is unique (c.f. [5, Theorem 3.2]). We denote its associated potential by Φ 0 and recall the relationship which solves the classical stationary system of self-gravitating particles (c.f. [3,15]), Now consider the energy functional associated to (1) where F is the Fisher information given by It is easy to see that E is weakly lower semicontinuous, strictly convex and coercive on P = n ∈ P | √ n ∈ H 1 (Ω) . Indeed, E 0 is equivalent to the functional which is uniformly bounded from below for all Φ ∈ H 1 0 (Ω) due to Moser [13]. Therefore, it attains a unique minimum n in the set P , and in particular in P. respectively. Then there exists (n * , Φ * , F * ) ∈ X × X × R, which solves the classical self-gravitation system (16) such that the following convergences hold for → 0 + : Furthermore, if n 0 is a unique minimizer of E 0 , then n * ≡ n 0 and Φ * ≡ Φ 0 . Proof. We begin by showing that { √ n } ⊂ H 1 (Ω) is bounded. Indeed, since n is a minimum of E , we have but E 0 (n 0 ) ≤ E 0 (n ) since n 0 is a minimum of E 0 . Hence, F(n ) ≤ F(n 0 ) for all > 0, which was to be shown. We can then extract a subsequence, denoted again by n , such that √ n √ n * in H 1 (Ω) and n → n * in L 2 (Ω) for some n * ∈ P, where the second convergence follows from the compact Sobolev embedding H 1 (Ω) → L 4 (Ω). Furthermore, we have which implies that E 0 (n 0 ) = lim →0+ E (n ). On the other hand, by the weakly lower L 2 (Ω)-semicontinuity of the functional E 0 , Therefore, n * is a minimizer of E 0 . The strong convergence Φ → Φ * := (−∆ 0 ) −1 n * in H 1 (Ω) follows easily from the strong convergence n → n * in L 2 (Ω), due to the continuity of the solution operator (−∆ 0 ) −1 . Now consider the Euler-Lagrange equation associated to E , i.e., the variational formulation of (1a), given by Similarly we have the Euler-Lagrange equation associated to E 0 , given by where F * ∈ R is the Lagrange multiplier for the constraint Ω n * dx = 1 with F * = − log e σΦ * L 1 (Ω) , since log n * = σΦ * + F * ∈ L ∞ (Ω).
Therefore, by testing the former variational formulation with ϕ = √ n and the latter with ϕ = n , and taking the difference of the resulting equations, we obtain where we have used the fact that F , F * ∈ R and Ω n dx = 1. Due to the convergences derived above, we conclude that F → F * in R. At this point, it is easy to see from (17) that the pair (u * , Φ * ) solves (16a). To see that it also solves (16b), we notice that Φ * is a minimizer of G and that (16b) is simply the Euler-Lagrange equation associated to G. Hence, (n * , Φ * ) ∈ X × X is indeed a solution of (16). Furthermore, from the representation u = log n −F we obtain the strong convergence u → σΦ * in L 2 (Ω). Indeed, by taking the difference of the two representations, multiplying the resulting equation with (u − σΦ * ) and integrating over Ω, we obtain which clearly yields the required convergence. Finally, if the minimizer n 0 is unique, then n * ≡ n 0 and consequently Φ * ≡ Φ 0 .
6. Numerical Simulations. In this section we show two numerical simulations that validate the theoretical results obtained above. In both cases, we considered the unit disk D ⊂ R 2 with scaled Planck constant = 1 × 10 −3 . An adaptive finite element method was used to solve the coupled problem (3) for (u, V ) iteratively. The outer iteration consist of a Picard iteration procedure for the potential V and the inner iteration consists of Newton's method to solve (3a) for the quasi potential u.
Writing the left-hand side of the inequality above in terms of ψ(G k (u)) and using the first inequality in (18) for half of the third term, we get Now decompose the right-hand side as follows For k 0 ≥ 4, we can absorb J 3 into the left-hand side of (19). As for J 1 , we have we can use the Hölder, interpolation and Young inequalities, along with the Sobolev embedding H 1 → L 2d d−2 to obtain So by choosing c 1 sufficiently small and k 0 sufficiently large, J 2 may also be absorbed into the left-hand side of (19), leaving us with for suitable constants δ i > 0, i = 1, 2. Here we used |s| ≤ |φ(s)| for |s| ≥ 0 and the fact that |G k (u)| = |u| − k a.e. on A k . Now if h > k > k 0 , then A h ⊂ A k and for arbitrary m ∈ N, From the definition of ζ we finally obtain the uniform bound |u| ≤ K 1 a.e. on Ω, which gives the bounds required.
A.2. Proof of Lemma 6.2. Following the proof of Lemma 6.1 we consider ϕ k = φ(G k (u)) for k ≥ k 0 (h) = max{1, 4h}, with a suitable h chosen later, as a test function in (6) to obtain This time we decompose the right-hand side as follows For each of the J i , i ∈ {1, 2, 3}, we have the bounds where we used the Sobolev embedding H 1 → L 2d d−2 and the second inequality in (18) for J 2 . By absorbing J 2 and J 3 into the left-hand side, we obtain δ ψ(G k (u)) 2 where, for h large enough, Hence ψ(G k (u)) ∈ H 1 (Ω) and ∇G k (u) ∈ [L 2 (Ω)] 2 for any k ≥ k 0 (h).