# A Two Species Thomas-Fermi Model for Stellar Ground States

Parker Hund1

###### Abstract

We present an energy functional for a Thomas-Fermi type two-fluid model of a
self-gravitating non-rotating charged body, with a non-relativistic kinetic
energy. We prove that, under certain conditions on the total number of
positively charged and negatively charged particles, a minimizer exists and
both fluids have compact support. We prove the same result for special
relativistic kinetic energy, assuming further conditions on the total number
of particles. In the non-relativistic kinetic energy case, we further prove
the uniqueness of the minimizer, as well as present results relating the
general shape of the minimizer to the total number of particles.

##  1 Introduction

The Thomas-Fermi model as an approximate model was created for the multi-
particle Schrödinger equation [1], [2]. In its molecular formulation, it is
concerned with minimizing the energy functional

$F(\rho)=C\int_{\mathbb{R}^{3}}\rho^{5/3}(x)d^{3}x+\int_{\mathbb{R}^{3}}\rho(x)V(x)d^{3}x+\frac{1}{2}q^{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho(x)\rho(y)}{|x-y|}d^{3}xd^{3}y$
(1)

under certain constraints on $\rho$. Here $\rho$ represents the electrons,
treated as fluid, which are in a potential field $V(x)$ created by fixed
atomic nuclei. The first term in the functional represents the kinetic energy
of the electrons while the last term represents the energy associated with the
electrical repulsion between electrons. A very detailed discussion of the
exact formulation of this model as well as a mathematically rigorous
discussion of its minimizers can be found in [11], [6].

How accurate this approximation to the energy of the ground state is was
unclear at the time, but it is now known that the Thomas-Fermi energy is the
leading order term, in the number of particles, of the ground state energy of
an atom [3]. As an approximate model it does have important limitations,
however. Most well known is the result of Teller that the theory does not
predict stable molecules. That is, one can lower the energy of a molecule by
breaking it apart into constituent atoms and moving them far apart from each
other [4]. This property was instrumental in Lieb and Thirring’s proof of
stability of non-relativistic matter by allowing for a reduction of “stability
of the second kind” to “stability of the first kind” [5]; that is, one reduces
the problems to atoms instead of molecules. A discussion of this idea and how
it appears in different proofs for the stability of matter can be found in
[7].

In this paper, we apply the Thomas-Fermi approximation to study the structure
of stellar objects instead of microscopic matter. This is somewhat similar to
the recent applications of Thomas-Fermi theory found in [12] and [13],
although those models are not presented in a mathematically rigorous way (nor
is their aim to be) and exhibit substantial differences when compared to ours.
It is more accurate to say that we follow in the footsteps of the effective
single fluid model for white dwarfs presented by Chandrasekhar [19] and the
results of Lieb and Yau in [17].

We essentially take Chandrasekhar’s model and forget the local neutrality
assumption. That is, he assumed that the negative charge of the electrons
would everywhere balance the positive charge of the various types of nuclei.
This allowed him to reduce the problem to finding just the electron density as
the density of the nuclei were assumed to be in constant proportion. We do not
make this assumption, but we do make the strong assumption that our stellar
object is made only of electrons and protons. This is necessary since the
Thomas-Fermi model applies to fermions but not bosons. Although a hydrogen
nucleus is a fermion, the other common nuclei like helium, carbon, or oxygen
found in stellar objects are all bosons. Therefore we must make the
restriction to only protons and electrons.

The result is a two species generalization of Thomas-Fermi model which treats
both the protons and electrons as fluids. This model can be used to describe
the ground state of a white dwarf star. We will, however, not dwell on the
physical motivation and justifications for this model: for a more detailed
introduction to its motivation, see [9] or [10].

The two fluid model of course has similarities to the single fluid Thomas-
Fermi model. But, because gravity will play an important role, it also shares
similarities to models of self-gravitating and rotating fluids. Discussions of
these closely related models can be found in [14], [15], [16], [17], and [18].
For reference, the functional used in [18] is given as

$F(\rho)=\int_{\mathbb{R}^{3}}A(\rho(x))d^{3}x-\int_{\mathbb{R}^{3}}\rho(x)J(x)d^{3}x-\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho(x)\rho(y)}{|x-y|}d^{3}xd^{3}y$
(2)

where $A(s)=s\int_{0}^{s}f(t)t^{-2}dt$, $f(\rho(x))=p(x)$ the equation of
state when $p(x)$ is pressure, and $J$ is the angular velocity.

The two fluid model is in some sense a combination of the single fluid Thomas-
Fermi model and the self-gravitating fluid models, something we will again
note in section 2, and in fact we will simply adapt some of the techniques
used for the rotating models in our study.

The outline of the paper is as follows. In section 2, we will present and
discuss the model. In section 3 we prove the existence of minimizers given
certain bounds on the ratio of the particles. In section 4 we present a
different proof under stricter conditions giving existence with compactly
supported minimizers. Section 5 presents an analysis of the Euler-Lagrange
equations of the energy functional. Section 6 relates the results of sections
3 and 4 to the results of section 5, and proves the uniqueness of the
minimizers. Finally, we conclude in section 7.

##  2 Energy Functionals

###  2.1 The non-relativistic model

Chandrasekhar’s white dwarf model assumes a special relativistic kinetic
energy. But before we generalize that model to two species, let us first
consider the non-relativistic case; it already presents enough mathematical
difficulties to make it interesting. The energy functional for the two species
non-relativistic model is

$\displaystyle E(\rho_{e},\rho_{p})$
$\displaystyle=\frac{3h^{2}}{40m_{e}}\left(\frac{3}{\pi}\right)^{2/3}\int_{\mathbb{R}^{3}}\rho^{5/3}_{e}(x)d^{3}x+\frac{3h^{2}}{40m_{p}}\left(\frac{3}{\pi}\right)^{2/3}\int_{\mathbb{R}^{3}}\rho^{5/3}_{p}(x)d^{3}x$
(3)
$\displaystyle+\frac{q^{2}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}(x)-\rho_{e}(x))(\rho_{p}(y)-\rho_{e}(y))}{|x-y|}d^{3}yd^{3}x$
$\displaystyle-\frac{G}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}(x)+m_{e}\rho_{e}(x))(m_{p}\rho_{p}(y)+m_{e}\rho_{e}(y))}{|x-y|}d^{3}yd^{3}x$

where $h$ is Planck’s constant, $m_{e}$ is the electron mass, $m_{p}$ is the
proton mass, $G$ is Newton’s gravitation constant, and $q$ is the elementary
charge. As a notational convenience, we will define
$k_{e}:=\frac{3h^{2}}{40m_{e}}\left(\frac{3}{\pi}\right)^{2/3}$ and
$k_{p}:=\frac{3h^{2}}{40m_{p}}\left(\frac{3}{\pi}\right)^{2/3}$.
Mathematically, the exact values of these constants will not make much
difference in what follows, and one could ignore any physical meaning
mentioned in these introductory sections. We will however retain the
suggestive notation and try to keep computed quantities in a form so that we
can more easily make analogies to well known physical quantities for those
readers familiar with, for example, the Chandrasekhar limiting mass. Although
we have included no numerical results in this paper, it should be noted that
for there are some computational issues with using the true physical values; a
detailed discussion of this issue as well as numerical results can be found in
[9].

As mentioned above, the first two terms on the right side of (3) are kinetic
energy terms. These are derived from a phase space argument combined with the
Pauli exclusion principle and nonrelativistic kinetic energies per particle,
$\frac{p^{2}}{2m}$. A careful explanation of these integrals can be found in
[11] or Chapter X of [19]. This is the term of the energy functional which
does not hold for bosons: note that the Pauli exclusion principle was used in
its justification. The third and fourth integrals are readily seen to be the
electric and gravitational energies.

It is also easy to see the similarities between this model and (2) and (1),
and why we might say that this model is, mathematically at least, in some way
a combination of those two models. From (1) we have taken the repulsive
potential and from (2) we have taken the attractive potential.

We wish to minimize $E$ subject to the following constraints:
$\rho_{e},\rho_{p}\geq 0$, $\rho_{e},\rho_{p}\in L^{1}\cap L^{5/3}$ (which is
sufficient to make the third and fourth integrals finite),
$\int\rho_{e}=N_{e}$, and $\int\rho_{p}=N_{p}$. We will sometimes refer to $E$
as $E_{(N_{e},N_{p})}$ when we want to make a point about particular values of
$N_{e}$ and $N_{p}$, but for the most part we will only use $E$. As mentioned
in the introduction, we are considering applications to stellar structures, so
we have in mind $N_{e},N_{p}$ on the order of $10^{56}$ (higher would make our
nonrelativistic energy assumption questionable), although it makes no
difference mathematically. Proofs that this functional is well defined and is
bounded below can be found in, for example, [18].

We can compute the Euler-Lagrange equations with $\int\rho_{e}=N_{e}$, and
$\int\rho_{p}=N_{p}$ constraints to be

$\frac{5}{3}k_{p}\rho_{p}^{2/3}(x)+q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}-\rho_{e})(y)}{|x-y|}d^{3}y-Gm_{p}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}+m_{e}\rho_{e})(y)}{|x-y|}d^{3}y=\lambda_{p}$
(4)

and

$\frac{5}{3}k_{e}\rho_{e}^{2/3}(x)-q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}-\rho_{e})(y)}{|x-y|}d^{3}y-Gm_{e}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}+m_{e}\rho_{e})(y)}{|x-y|}d^{3}y=\lambda_{e}$
(5)

These will hold where $\rho_{p}>0$ and $\rho_{e}>0$ respectively. If the
densities were sufficiently regular, we could apply $-\Delta$ to (4) and (5)
to get

$-\frac{5}{3}k_{p}\Delta\rho^{2/3}_{p}+4\pi q^{2}(\rho_{p}-\rho_{e})-4\pi
Gm_{p}(m_{p}\rho_{p}+m_{e}\rho_{e})=0$ (6)

and

$-\frac{5}{3}k_{e}\Delta\rho^{2/3}_{e}-4\pi q^{2}(\rho_{p}-\rho_{e})-4\pi
Gm_{e}(m_{p}\rho_{p}+m_{e}\rho_{e})=0$ (7)

again valid where $\rho_{p}>0$ and $\rho_{e}>0$. Note that these are the
equations of hydrostatic equilibrium for a self-gravitating, charged two fluid
model if one use a polytropic equation of state of index 3/2. Since both
functions are positive on their support, we define
$\nu_{e,p}=\rho^{2/3}_{e,p}$. Then we rewrite (6) and (7) as

$\frac{5}{3}k_{p}\Delta\nu_{p}=4\pi q^{2}(\nu^{3/2}_{p}-\nu^{3/2}_{e})-4\pi
Gm_{p}(m_{p}\nu^{3/2}_{p}+m_{e}\nu^{3/2}_{e})$ (8)

and

$\frac{5}{3}k_{e}\Delta\nu_{e}=-4\pi q^{2}(\nu^{3/2}_{p}-\nu^{3/2}_{e})-4\pi
Gm_{e}(m_{p}\nu^{3/2}_{p}+m_{e}\nu^{3/2}_{e})$ (9)

So, aside from the question of regularity, we have reduced finding the
critical points of (3) to finding solutions to this nonlinear elliptic system.

Positive solutions to systems of this form have been studied extensively; see
[20] for an overview. The common techniques to find solutions are usually
extensions of the techniques applied to find positive solutions of the scalar
versions. For example, [21] generalizes the results of [23], [24] generalizes
[25] and [26], and [22] generalizes [27]. These techniques cannot, at least
straightforwardly, be applied to our system. One of the main issues is that
many of the techniques use a fixed point iteration on a function space of
positive functions. In our system, we cannot force the right hand sides of (8)
and (9) to be positive, preventing us from applying the maximum principle in a
crucial way.

But the more serious issue is that these techniques are simply not addressing
the question we are interested in. These papers are overwhelmingly concerned
with finding strictly positive solutions with Dirichlet boundary conditions on
an arbitrary domain. While for spherical domains there do exist solutions
satisfying these criterion, in fact easily describable as we will show below,
we are also interested in solutions where one density has support strictly
contained in the support of the other density. As we will show, there are only
special pairs $(N_{e},N_{p})$ such that the densities have the same support.

###  2.2 The Special Relativistic Model

Now we may turn to the direct generalization of the Chandrasekhar model. As
mentioned above, the kinetic energy terms $k_{e}\int\rho_{e}^{5/3}$ and
$k_{p}\int\rho_{p}^{5/3}$ are derived by assuming a nonrelativistic kinetic
energy per particle of $\frac{p^{2}}{2m}$. If instead one uses the special-
relativistic kinetic energy of $mc^{2}((1+p^{2}/2m)^{1/2}-1)$, one gets
kinetic energy terms of the form

$\frac{\pi
m_{f}^{4}c^{5}}{3h^{3}}\int_{\mathbb{R}^{3}}A\left(\frac{h}{m_{f}c}\left(\frac{3}{8\pi}\right)^{1/3}\rho_{f}^{1/3}(x)\right)d^{3}x$
(10)

where

$A(z)=8z^{3}[(z^{2}+1)^{1/2}-1]-z(2z^{2}-3)(z^{2}+1)^{1/2}-3\sinh^{-1}(z)$
(11)

and $f=e$ or $p$, see [19] Chapter X for the derivation of exactly this
expression. This is however not always the best form of the kinetic energy to
work with. Sometimes a better form is

$m_{f}c^{2}\int_{\mathbb{R}^{3}}\int_{0}^{\rho_{f}(x)}\sqrt{1+\left(\frac{3}{\pi}\right)^{2/3}\left(\frac{h}{2m_{f}c}\right)^{2}\theta^{2/3}}d\theta
d^{3}x.$ (12)

So to express the special relativistic energy functional, we only need to
replace the first two terms in (3) with their analogues from (12). We will
call this energy functional $E^{S}$. Then we make the same assumption of
spherical symmetry and seek to minimize $E^{S}$ over the same set as before.
If we compute the Euler-Lagrange equations of $E^{S}$ and assume high enough
regularity to apply $-\Delta$, we find as analogues to (6) and (7)

$-m_{p}c^{2}\Delta\sqrt{1+\left(\frac{3}{\pi}\right)^{2/3}\left(\frac{h}{2m_{p}c}\right)^{2}\rho_{p}^{2/3}}+4\pi
q^{2}(\rho_{p}-\rho_{e})-4\pi Gm_{p}(m_{p}\rho_{p}+m_{e}\rho_{e})=0$ (13)

and

$-m_{e}c^{2}\Delta\sqrt{1+\left(\frac{3}{\pi}\right)^{2/3}\left(\frac{h}{2m_{e}c}\right)^{2}\rho_{e}^{2/3}}-4\pi
q^{2}(\rho_{p}-\rho_{e})-4\pi Gm_{e}(m_{p}\rho_{p}+m_{e}\rho_{e})=0.$ (14)

If, following [19] Chapter XI, we set
$y_{f}^{2}=1+\left(\frac{3}{\pi}\right)^{2/3}\left(\frac{h}{2m_{f}c}\right)^{2}\rho_{f}^{2/3}$
where $f$ is $p$ or $e$, these equations can be rewritten as

$\Delta
y_{p}=\frac{4\pi}{m_{p}c^{2}}(q^{2}-Gm_{p}^{2})\left(\frac{\pi}{3}\right)\left(\frac{2m_{p}c}{h}\right)^{3}(y_{p}^{2}-1)^{3/2}-\frac{4\pi}{m_{p}c^{2}}(q^{2}+Gm_{p}m_{e})\left(\frac{\pi}{3}\right)\left(\frac{2m_{e}c}{h}\right)^{3}(y_{q}^{2}-1)^{3/2}$
(15)

and

$\Delta
y_{e}=-\frac{4\pi}{m_{e}c^{2}}(q^{2}+Gm_{p}m_{e})\left(\frac{\pi}{3}\right)\left(\frac{2m_{p}c}{h}\right)^{3}(y_{p}^{2}-1)^{3/2}+\frac{4\pi}{m_{e}c^{2}}(q^{2}-Gm_{q}^{2})\left(\frac{\pi}{3}\right)\left(\frac{2m_{e}c}{h}\right)^{3}(y_{q}^{2}-1)^{3/2},$
(16)

these equations valid where $\rho_{p}>0$ and $\rho_{e}>0$, respectively. These
are clearly much easier to work with than (13) and (14). In his single species
model, Chandrasekhar is able to go a step further and introduce a change of
variables to scale out the constants and essentially normalize the equation.
He is left with an equation (assuming spherical symmetry) of the form

$\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{dy}{dr}\right)=-\left(y^{2}-\frac{1}{y_{0}^{2}}\right)^{3/2}$
(17)

where $y_{0}$ is a function of the central density. Rewriting the equation
this way allows him to write certain quantities, such as the total mass of the
configuration, with a simple expression. This is what allows him to write his
well known limiting mass of possible configurations in a closed form.

Since we are working with a system, there is no scaling which will eliminate
all the constants of (15) and (16) to get a system analogous to (17). We
therefore must work with (15) and (16) and will not focus on the question of
limiting mass, although a few of our results will shed some light on it.

###  2.3 Zero-gravity limit

Before we turn to proving the existence of minimizers for (3) and its special
relativistic form, it should be noted that, although the gravitational force
is very small in comparison to the electric force, setting $G=0$ trivializes
the problem, at least in the non-relativistic case. Define

$\displaystyle E^{0}(\rho_{e},\rho_{p})$
$\displaystyle=k_{e}\int_{\mathbb{R}^{3}}\rho^{5/3}_{e}(x)d^{3}x+k_{p}\int_{\mathbb{R}^{3}}\rho^{5/3}_{p}(x)d^{3}x$
(18)
$\displaystyle+\frac{q^{2}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}(x)-\rho_{e}(x))(\rho_{p}(y)-\rho_{e}(y))}{|x-y|}d^{3}yd^{3}x$

We show that there is only a single local critical point, which happens to be
the global minimum. The first two terms of $E^{0}$ are clearly positive, and
the proof that the last term is positive can be found in Theorem II.6 in [11].
Therefore, $E^{0}$ has a global minimum of zero, obtained only when $\rho_{e}$
and $\rho_{p}$ are identically zero. But beyond this, for any fixed $N_{p}$
and $N_{e}$, $E^{0}$ has an infimum of zero. To see this, consider any
$\rho_{p}$ and $\rho_{e}$ in the set of admissible functions. Then let
$\rho^{\lambda}_{e}(s)=\rho_{e}(s/\lambda)/\lambda^{3}$ and
$\rho^{\lambda}_{p}(s)=\rho_{p}(s/\lambda)/\lambda^{3}$. We have

$\int_{\mathbb{R}^{3}}\rho_{p}^{\lambda}(s)d^{3}s=\int_{\mathbb{R}^{3}}\frac{\rho_{p}(s/\lambda)}{\lambda^{3}}d^{3}s=\int_{\mathbb{R}^{3}}\rho_{p}(s)d^{3}s=N_{p}$
(19)

so $\rho_{p}^{\lambda}$ is admissible for $0<\lambda<\infty$. We also have

$\int_{\mathbb{R}^{3}}(\rho_{p}^{\lambda})^{5/3}(s)d^{3}s=\int_{\mathbb{R}^{3}}\frac{\rho_{p}^{5/3}(s/\lambda)}{\lambda^{5}}d^{3}s=\int_{\mathbb{R}^{3}}\frac{\rho_{p}^{5/3}(s)}{\lambda^{2}}d^{3}s\rightarrow
0\text{ as }\lambda\rightarrow\infty$ (20)

These statements are also true for $\rho_{e}$. We also have for the third term
in $E^{0}$:

$\displaystyle\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho^{\lambda}_{p}(s)-\rho^{\lambda}_{t}(s))(\rho^{\lambda}_{p}(t)-\rho^{\lambda}_{e}(t))}{|s-t|}d^{3}sd^{3}t$
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}(s/\lambda)-\rho_{t}(s\lambda))(\rho_{p}(t/\lambda)-\rho_{e}(t/\lambda))}{|s-t|\lambda^{6}}d^{3}sd^{3}t$
(21)
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}(s)-\rho_{t}(s))(\rho_{p}(t)-\rho_{e}(t))}{|s-t|\lambda}d^{3}sd^{3}t$
(22) $\displaystyle\rightarrow 0\text{ as }\lambda\rightarrow\infty$

This process is physically analogous to spreading any given density out. We
can then conclude that the infimum of $E^{0}$ for any fixed pair
$(N_{e},N_{p})$ is zero, and not obtained.

Beyond this, $E^{0}$ is also jointly convex, so the global minimum is the only
critical point. To see this note that the first two terms of $E^{0}$ are
strictly convex on $[0,\infty]$ since $x^{5/3}$ is. For the last term,
consider

$\hat{E}^{0}(\rho_{p},\rho_{e})=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}(s)-\rho_{t}(s))(\rho_{p}(t)-\rho_{e}(t))}{|s-t|}d^{3}sd^{3}t$
(23)

Let $u_{e}$ and $u_{p}$ be admissible functions and define
$\hat{E}^{0}(\alpha)=\hat{E}_{G}(\alpha\rho_{p}+(1-\alpha)u_{p},\alpha\rho_{e}+(1-\alpha)u_{e})$.
Then we can compute

$\frac{d^{2}\hat{E}^{0}}{d\alpha^{2}}|_{\alpha=0}=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(u_{p}(s)-u_{t}(s))(u_{p}(t)-u_{e}(t))}{|s-t|}d^{3}sd^{3}t$
(24)

As was stated above, this quantity is positive and only zero if $u_{e}=u_{p}$,
so $E^{0}$ is jointly convex. So one might have considered using a
perturbative expansion in $G$ around $G=0$ to approach finding solutions. But
since $E^{0}$ has only trivial critical points this does not simplify the
analysis; see [30] for this approach.

##  3 Minimizers of the Energy Functional

This section generally follows the outline of [31].

######  Theorem 3.1.

If

$\frac{1-\frac{Gm_{p}^{2}}{q^{2}}}{1+\frac{Gm_{e}m_{p}}{q^{2}}}\leq\frac{N_{e}}{N_{p}}\leq\frac{1+\frac{Gm_{e}m_{p}}{q^{2}}}{1-\frac{Gm_{q}^{2}}{q^{2}}},$
(25)

there is a pair $(\rho_{e},\rho_{p})$ which minimizes $E$ such that
$\rho_{e},\rho_{p}\geq 0$, $\rho_{e},\rho_{p}\in L^{1}\cap L^{5/3}$,
$\int\rho_{e}=N_{e}$, $\int\rho_{p}=N_{p}$, and $\rho_{e},\rho_{p}$ are both
radially symmetric.

Define $W^{N_{e},N_{p}}$ to be the set of pairs of functions
$(\rho_{e},\rho_{p})$ such that $\rho_{f}\geq 0$, $\rho_{f}\in L^{1}\cap
L^{5/3}$, $\rho_{f}$ are radially symmetric, and $\int\rho_{f}\leq N_{f}$.
Also, let $\mathcal{E}^{N_{e},N_{p}}$ be the minimum energy of $E$ over
$W^{N_{e},N_{p}}$; this could possibly be $-\infty$. To avoid the clunky
notation, we will always fix $N_{e}$ and $N_{p}$ and just write $W$ and
$\mathcal{E}$.

Sometimes we will also use the notation

$B\rho(x):=\int_{\mathbb{R}^{3}}\frac{\rho(y)}{|x-y|}d^{3}y.$ (26)

First, we want to show

######  Lemma 3.2.

For any $(N_{e},N_{p})$ satisfying $(\ref{5/3bounds})$, $\mathcal{E}<0$.

###### Proof.

It is easy to see that this is true in the case that $N_{e}=N_{p}$. For if we
take any configuration in $W$ such that $\rho_{e}=\rho_{p}$, the electric
energy is zero and we are left with

$E(\rho_{e},\rho_{p})=(k_{e}+k_{p})\int_{\mathbb{R}^{3}}\rho^{5/3}_{e}(x)d^{3}x-\frac{G}{2}(m_{e}+m_{p})^{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho_{e}(x)\rho_{e}(y)}{|x-y|}d^{3}yd^{3}x$
(27)

Using the same “spreading out” argument we used when we considered the zero-
gravity case, we see that as $\lambda\rightarrow\infty$, the kinetic energy
terms go to zero like $1/\lambda^{2}$ while the gravitational energy goes to
zero like $1/\lambda$. Then for large enough $\lambda$,
$E(\rho_{e}^{\lambda},\rho_{p}^{\lambda})<0$.

Therefore, because each of the terms in $E$ is in some sense continuous in
$(\rho_{e},\rho_{p})$ (without trying to make this precise in any way, we just
use the intuition), there is some interval of values for $N_{e}/N_{p}$
containing 1 such that $\mathcal{E}<0$. We determine the size of this interval
by finding example functions in $W$.

Suppose that for a fixed $(N_{e},N_{p})$ such that the ratio $N_{e}/N_{p}$
satisfies $(\ref{5/3bounds})$, we have a pair $(\rho_{e},\rho_{p})$ both with
compact support such that $E(\rho_{e},\rho_{p})<0$. We claim we can add a
small positive or negative charge outside of this configuration and the energy
will not increase. Formally, the idea is to see how many “test particles” we
can bring in from infinity and still have negative energy. Of course, this is
not a rigorous argument, so we now make this idea rigorous.

Let us assume we are adding positive charge. Since our space $W$ consists of
spherically symmetric functions, we must add our charge in a spherically
symmetric configuration. We may assume the configuration to which we are
adding has compact support, suppose it is contained in a ball of radius
$R_{0}$. Then add to $(\rho_{e},\rho_{p})$ a particle density described by the
function $h=\epsilon\chi_{S}$ for $S=\\{x|R_{0}<|x|<R_{0}+\eta\\}$. Then we
compute

$\displaystyle E(\rho_{e},\rho_{p}+g)-E(\rho_{e},\rho_{p})$
$\displaystyle=k_{p}\epsilon^{5/3}\frac{4}{3}\pi\left[(R_{0}+\eta)^{3}-R_{0}^{3}\right]+\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}\int_{S}\int_{S}\frac{1}{|x-y|}d^{3}yd^{3}x$
$\displaystyle+\epsilon\int_{S}\int_{\mathbb{R}^{3}}\frac{q^{2}(\rho_{p}-\rho_{e})-Gm_{p}(m_{p}\rho_{p}+m_{e}\rho_{e})}{|x-y|}$
(28)
$\displaystyle=k_{p}\epsilon^{5/3}\frac{4}{3}\pi\left[3R_{0}^{2}\eta+3R_{0}\eta^{2}+\eta^{3}\right]$
$\displaystyle+\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}4\pi\int_{R_{0}}^{R_{0}+\eta}r^{2}\left[\frac{4}{3}\frac{\pi}{r}\left[r^{3}-R_{0}^{3}\right]+2\pi\left[(R_{0}+\eta)^{2}-r^{2}\right]\right]dr$
$\displaystyle+4\pi\epsilon\int_{R_{0}}^{R_{0}+\eta}\left[q^{2}(N_{p}-N_{e})-Gm_{p}(m_{p}N_{p}+m_{e}N_{e})\right]rdr$
(29)
$\displaystyle=k_{p}\epsilon^{5/3}\frac{4}{3}\pi\left[3R_{0}^{2}\eta+3R_{0}\eta^{2}+\eta^{3}\right]$
$\displaystyle+\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}4\pi\left[4\pi
R_{0}^{3}\eta^{2}+\frac{16\pi}{3}R_{0}^{2}\eta^{3}+\frac{8\pi}{3}R_{0}\eta^{4}+\frac{3\pi}{5}\eta^{5}\right]$
$\displaystyle-4\pi\epsilon
N_{e}(q^{2}-Gm_{p}^{2})\left(B_{p}-\frac{N_{p}}{N_{e}}\right)(2R_{0}\eta+\eta^{2})$
(30)

where $B_{p}$ is the multiplicative inverse of the left side of (25). If we
assume $\left(B_{p}-\frac{N_{p}}{N_{e}}\right)>0$, the last term is strictly
negative. Therefore, we may take $\epsilon>0$ small enough so that the
difference is negative. So we can add some positive amount of charge to the
configuration and the energy will become more negative. The same will be true
of adding a small negative amount of charge using the strict version of the
right half of (25).

We intend to treat this as an iterative procedure, and note that if
$N_{e}=N_{p}$, we can assume that the densities have compact support; the size
of the support makes no difference for the argument that
$\mathcal{E}^{N_{e},N_{e}}<0$. The issue remaining is whether we can saturate
(25) using this procedure.

To show that we can, we let $\eta=\epsilon$. Then if we have only added a
finite amount of charge to the configuration, we can also be assured the
resulting configuration has finite radius. So no matter how much charge we
add, we can always apply the procedure again as long as the strict version of
(25) is satisfied. This will give us (25). So we need to show that the
procedure can be carried out if $\eta=\epsilon$, which is to say that we can
find an $\epsilon$ so that in this case (4) is negative.

The leading order in $\epsilon$ parts of the three terms of (4) behave as
$\epsilon^{8/3}R_{0}^{2}$, $\epsilon^{3}R_{0}^{4}$, and $\epsilon^{2}R_{0}$,
respectively. So if we take $\epsilon$ to be $o(R_{0}^{-4})$, we can make the
first two terms smaller in absolute value than the last, negative term. ∎

Let

$u_{f}^{n}(x)=\int_{\mathbb{R}^{3}}\frac{\rho_{f}^{n}(y)}{|x-y|}d^{3}y$ (31)

for $f=e,p$. Then due to the uniform bound on the $L^{4/3}$ norm, we have that
$\nabla u_{f}^{n}$ is bounded in $L^{q}$ for $3/2<q\leq 12/5$ and $u_{f}^{n}$
is bounded in $L^{p}$ for $3<p\leq 12$, see [14] proposition 6. Note also that
$u_{f}^{n}$ is radially symmetric. Since $\rho_{f}^{n}$ is uniformly bounded
in the $L^{4/3}$ norm, we can apply Banach-Alaoglu and find a weakly
convergent subsequence. From this subsequence, extract a further subsequence
converging weakly in $L^{6/5}$. Let $(\rho_{e},\rho_{p})$ be the functions to
which $(\rho_{e}^{n},\rho_{p}^{n})$ converge. We can extract yet a further
subsequence so that $(u_{e}^{n},u_{p}^{n})$ converge weakly in $L^{6}$ and
a.e. to $(u_{e},u_{p})$.

The goal is to prove that $(\rho_{e},\rho_{p})$ is a minimizer of $E$ over $W$
and $\int\rho_{p}=N_{p}$ and $\int\rho_{e}=N_{e}$. To see the former, we show
that $E$ is weakly lower semicontinuous. First note that

$H(g):=\int g^{5/3}\leq\liminf\int g_{n}^{5/3}$ (32)

for the same reasons as given in [31] (in that paper, instead of $g^{5/3}$
there is given a general $j(\rho)$ nonnegative, continuous, and convex). So we
only need to focus on proving that the two other terms in the energy are lower
semicontinuous. This result follows by proving that both $(u_{p}^{n})$ and
$(u_{e}^{n})$ are relatively compact in the $(L^{6/5})^{*}=L^{6}$ topology and
so some subsequence converges to $(u_{e},u_{p})$. Coupled with the weak
$L^{6/5}$ convergence of $\rho_{f}^{n}$, we then have

$\int_{\mathbb{R}^{3}}(m_{p}\rho_{p}^{n}(x)+m_{e}\rho_{e}^{n}(x))(m_{p}u_{p}^{n}(x)+m_{e}u_{e}^{n}(x))d^{3}x\rightarrow\int_{\mathbb{R}^{3}}(m_{p}\rho_{p}(x)+m_{e}\rho_{e}(x))(m_{p}u_{p}(x)+m_{e}u_{e}(x))d^{3}x,$
(33)

and

$\int_{\mathbb{R}^{3}}(\rho_{p}^{n}(x)-\rho_{e}^{n}(x))(u_{p}^{n}(x)-u_{e}^{n}(x))d^{3}x\rightarrow\int_{\mathbb{R}^{3}}(\rho_{p}(x)-\rho_{e}(x))(u_{p}(x)-u_{e}(x))d^{3}x.$
(34)

To see that $(u_{f}^{n})$ is relatively compact in $L^{6}$, we first apply
Proposition II.1 from [31] to obtain the bound

$|u_{f}^{n}(x)|\leq C(N,p,q)\left(||\nabla
u_{f}^{n}||_{L^{p}}^{p^{\prime}/q+p^{\prime}}||u_{f}^{n}||_{L^{q}}^{q/q+p^{\prime}}\right)|x|^{-2p^{\prime}/(q+p^{\prime})}$
(35)

where $p^{\prime}$ is the Hölder conjugate of of $p$ and $C$ is a positive
constant. This allows us to conclude that $u_{f}^{n}\rightarrow 0$ as
$|x|\rightarrow\infty$ uniformly in $x$ and $n$. This is requirement v) of
Compactness Lemma 2 from [33]. Since we already have a.e. convergence of
$u_{f}^{n}$, it is straightforward to check that with $P^{n}(s)=s^{6}$ and
$Q^{n}(s)=s^{7}+\chi_{B_{1}}(s)s^{5}$ where $\chi$ is an indicator function,
the other four requirements are also satisfied. We can then apply the lemma to
conclude that

$\int_{\mathbb{R}^{3}}||u_{f}^{n}(x)|^{6}-|u_{f}(x)|^{6}|d^{3}x\rightarrow 0$
(36)

This gives strong convergence of $(u_{f}^{n})$ to $u_{f}$ in $L^{6}$.

Using (32), (33), and (34), we can conclude that $E$ is weakly lower
semicontinuous, and therefore

$E(\rho_{e},\rho_{p})\leq\liminf E(\rho_{e}^{n},\rho_{p}^{n})$ (37)

So we have a minimizer $(\rho_{e},\rho_{p})$. However, we can only so far
conclude that that $\int\rho_{p}\leq N_{p}$ and $\int\rho_{e}\leq N_{e}$; we
need to prove these are actual equalities.

We will use the following well known lemma, which can be found in, for
example, [14] as proposition 5:

######  Lemma 3.3.

Suppose $\rho\in L^{1}\cap L^{p}$. If $1<p\leq 3/2$, then $B\rho$ is in
$L^{r}$ for $3<r<3p/(3-2p)$, and

$||B\rho||_{r}\leq
c_{0}(||\rho||_{1}^{b}||\rho||_{p}^{1-b}+||\rho||_{1}^{c}||\rho||_{p}^{1-c})$
(38)

where $0<b<1$ and $0<c<1$.

######  Lemma 3.4.

The minimizer $(\rho_{e},\rho_{p})$ satisfies $\int\rho_{p}=N_{p}$ and
$\int\rho_{e}=N_{e}$.

###### Proof.

First, assume that both $\int\rho_{e}<N_{e}$ and $\int\rho_{p}<N_{p}$. We can
argue in a way similar to [31]. Let
$\rho_{f}^{\sigma}(x)=\rho_{f}(\frac{x}{\sigma^{1/3}})$. Then

$\int_{\mathbb{R}^{3}}\rho_{f}^{\sigma}(x)d^{3}x=\sigma\int_{\mathbb{R}^{3}}\rho(x)d^{3}x$
(39)

and

$\displaystyle E(\rho_{e}^{\sigma},\rho_{p}^{\sigma})$
$\displaystyle=\sigma\left(k_{e}\int_{\mathbb{R}^{3}}\rho^{5/3}_{e}d^{3}x+k_{p}\int_{\mathbb{R}^{3}}\rho^{5/3}_{p}d^{3}x\right)$
(40)
$\displaystyle+\frac{q^{2}\sigma^{5/3}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}(x)-\rho_{e}(x))(\rho_{p}(y)-\rho_{e}(y))}{|x-y|}d^{3}yd^{3}x$
$\displaystyle-\frac{G\sigma^{5/3}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}(x)+m_{e}\rho_{e}(x))(m_{p}\rho_{p}(y)+m_{e}\rho_{e}(y))}{|x-y|}d^{3}yd^{3}x$

Then we can compute

$\frac{dE}{d\sigma}(\rho_{e}^{\sigma},\rho_{p}^{\sigma})|_{\sigma=1}=K(\rho_{e},\rho_{p})+\frac{5}{3}V(\rho_{e},\rho_{p})=E(\rho_{e},\rho_{p})+\frac{2}{3}V(\rho_{e},\rho_{p})$
(41)

where we have called the kinetic energy terms $K$ and the potential energy
terms $V$, so that
$E(\rho_{e},\rho_{p})=K(\rho_{e},\rho_{p})+V(\rho_{e},\rho_{p})$. But we know
that $E(\rho_{e},\rho_{p})<0$, and since $K(\rho_{e},\rho_{p})>0$, we must
then have $V(\rho_{e},\rho_{p})<0$. Therefore,

$\frac{dE}{d\sigma}(\rho_{e}^{\sigma},\rho_{p}^{\sigma})|_{\sigma=1}<0,$ (42)

so we can increase $\sigma$ slightly so that we still have $\int\rho_{e}\leq
N_{e}$ and $\int\rho_{p}\leq N_{p}$ but the energy is decreased. This
contradicts $(\rho_{e},\rho_{p})$ being a minimizer.

This leaves us with the two cases such $\int\rho_{e}=N_{e}$ while
$\int\rho_{p}<N_{p}$ and $\int\rho_{e}<N_{e}$ while $\int\rho_{p}=N_{p}$.
These cases can be treated in the same way, so let us consider the former. If
both densities have bounded support, then we are in the same situation as we
considered in lemma (3.2): adding a small sphere of positive charge outside of
both supports will actually decrease the energy of the configuration. Then we
contradict that $(\rho_{e},\rho_{p})$ is a minimizer.

If either or both supports are unbounded, define $N_{f}(R)$ to be the number
of particles of type $f$ within a sphere of radius $R$. Then the function
$F(R)=q^{2}(N_{p}(R)-N_{e}(R))-Gm_{p}(m_{p}N_{p}(R)+m_{e}N_{e}(R))$ is
continuous, and approaches some negative value as $R\rightarrow\infty$. Let us
choose $R_{0}$ so that $F(R)<\beta<0$ for $R>R_{0}$. Now we add a spherical
shell of positive charge of width $\eta$, letting $h$ be as in lemma 3.2, and
compute the change in energy. We find that

$\displaystyle
E(\rho_{e},\rho_{p}+h)-E(\rho_{e},\rho_{p})=k_{p}\int_{S}((\rho_{p}+\epsilon)^{5/3}-\rho_{p}^{5/3})d^{3}x$
(43) $\displaystyle+$
$\displaystyle\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}4\pi\left[4\pi
R_{0}^{3}\eta^{2}+\frac{16\pi}{3}R_{0}^{2}\eta^{3}+\frac{8\pi}{3}R_{0}\eta^{4}+\frac{3\pi}{5}\eta^{5}\right]$
$\displaystyle+$ $\displaystyle
4\pi\epsilon\int_{R_{0}}^{R_{0}+\eta}\left[q^{2}(N_{p}(r)-N_{e}(r))-Gm_{p}(m_{p}N_{p}(r)+m_{e}N_{e}(r))\right]rdr$
$\displaystyle+$ $\displaystyle
4\pi\epsilon\int_{R_{0}}^{R_{0}+\eta}r^{2}\int_{|y|>R_{0}+r}\frac{q^{2}(\rho_{p}-\rho_{e})(y)-m_{p}G(m_{p}\rho_{p}+m_{e}\rho_{e})(y)}{|y|}d^{3}ydr$

The two differences from the situation we considered before are that now we
must be concerned with the kinetic energy term, and the final term,
corresponding to the interaction between the sphere of added particles and the
configuration outside this sphere. We want to show that the third term on the
right side is larger in absolute value than the other terms. By our assumption
on the size of $R_{0}$, we can say that this term is less than
$-4\pi\epsilon\eta R_{0}\beta$

To leading order in $\epsilon$, the kinetic energy term is $\pi
k_{p}\frac{20}{3}\epsilon\int_{R_{0}}^{R_{0}+\eta}\rho_{p}^{2/3}(r)dr$.
Because $\int\rho_{p}=N_{p}$, we must be able to find a set of total width in
the radial direction of $\eta$ and with radial distance at least $R_{0}$ such
that $\rho_{p}(r)\leq\frac{1}{r^{3}}$. We will assume for simplicity that this
set is $S$, the spherical shell, but it may be that $S$ needs to be broken up
into many unconnected shells. Then the kinetic energy term is less than $\pi
k_{p}\frac{20}{3}\epsilon\frac{\eta}{R_{0}}$, clearly smaller than the third
term for large enough $R_{0}$.

The second term is smaller than the third by taking both $\epsilon$ and $\eta$
small. For the last term, since both $\rho_{p}$ and $\rho_{e}$ are integrable,
one can choose $R_{0}$ large enough so that
$\int_{|y|>R_{0}}q^{2}(\rho_{p}+\rho_{e})(y)<\hat{\epsilon}$ for any chosen
$\hat{\epsilon}$. Then we can say that the last term is bounded by
$\frac{4\pi\epsilon}{3R_{0}}(3R_{0}\eta^{2}+3\eta
R^{2}_{0}+\eta^{3})\hat{\epsilon}$. The leading order term in $\epsilon$ is
then $4\pi\epsilon\eta R_{0}\hat{\epsilon}$. Compare this to the leading order
term for the negative third term, $4\pi\epsilon^{2}R_{0}\beta$. Then we make
$R_{0}$ at least large enough so that $\hat{\epsilon}<|\beta|$, and then
choose $\epsilon$ so that the negative third term dominates the positive
terms.

This contradicts that we have found a minimizer, and we find that we have
$\int\rho_{e}=N_{e}$ and $\int\rho_{p}=N_{p}$. Note that we have assumed
continuity of the minimizer for simplicity. This will independently be proven
below. ∎

This completes the proof of the theorem.

###  3.1 The Special Relativistic Case

The same proof also works for the special relativistic case with a few minor
changes which we will discuss below. However, recall that in the special
relativistic case we do not know that the solutions are decreasing, so our
result is not as strong. Here we can redefine $W$ so that it contains
functions in $L^{1}\cap L^{4/3}$ instead of $L^{1}\cap L^{5/3}$. We then have

######  Theorem 3.5.

If

$\frac{1-\frac{Gm_{p}^{2}}{q^{2}}}{1+\frac{Gm_{e}m_{p}}{q^{2}}}\leq\frac{N_{e}}{N_{p}}\leq\frac{1+\frac{Gm_{e}m_{p}}{q^{2}}}{1-\frac{Gm_{q}^{2}}{q^{2}}},$
(44)

and $[m_{p}N_{p}+m_{e}N_{e}]^{2/3}<\frac{\pi
2^{2/3}hc\left(\frac{3}{8\pi}\right)^{4/3}}{GKm_{p}^{4/3}}$, then there is a
pair $(\rho_{e},\rho_{p})$ which minimizes $E^{s}$ in $W$.

The constant $K$ here is the same constant which appears in [31] equation
(5’). To save space, let $\alpha_{f}$ denote
$\frac{h}{m_{f}c}\left(\frac{3}{8\pi}\right)^{1/3}$.

  1. 1.

In the proof that the energy of the minimizer is negative if (44) is
satisfied, the first step was to consider the $N_{p}=N_{e}$ case. The
“spreading out” argument does not work as nicely since $g(z)$ is not a
homogeneous polynomial. But if we take a $\rho_{f}$ which is bounded and again
say that $\rho_{f}^{\lambda}(x)=\rho_{f}(x/\lambda)/\lambda^{3}$, we can say
that, since $g(z)$ behaves like $z^{5}$ as $z\rightarrow 0$, we have

$\displaystyle\int_{\mathbb{R}^{3}}g(\alpha_{f}(\rho_{f}^{\lambda}(x))^{1/3})d^{3}x$
$\displaystyle=\int_{\mathbb{R}^{3}}g\left(\frac{\alpha_{f}}{\lambda}\left(\rho_{f}\left(\frac{x}{\lambda}\right)\right)^{1/3}\right)d^{3}x$
(45)
$\displaystyle\leq\int_{\mathbb{R}^{3}}M\left(\frac{\alpha_{f}}{\lambda}\right)^{9/2}\rho_{f}\left(\frac{x}{\lambda}\right)^{3/2}d^{3}x$
(46)
$\displaystyle=\frac{1}{\lambda^{3/2}}\int_{\mathbb{R}^{3}}M\alpha_{f}^{9/2}\rho_{f}(x)^{3/2}d^{3}x$
(47)

for large $\lambda$ and some $M$. It therefore decays in $\lambda$ faster than
the gravitational energy, and we can find a negative energy state.

We can use the same idea with the second part of the argument: since the
argument is that we can shrink the added density until the gravitational
energy dominates, the kinetic energy term will again behave like
$\epsilon^{5/3}$ for small $\epsilon$.

  2. 2.

The process of obtaining the uniform $L^{4/3}$ bounds on $\rho^{n}_{f}$ must
be changed to add the extra condition from the statement of theorem 3.5.
Before, we took the bound on the gravitational energy,

$C_{1}=2^{1/3}GKm_{p}^{4/3}(m_{p}N_{p}+m_{e}N_{e})^{2/3},$ (48)

and found a lower bound on $\int\rho^{5/3}$ involving
$(C_{1}+\epsilon)\int\rho^{4/3}$. This was to bound the negative gravitational
energy in terms of the positive kinetic energy. We were able to do this
regardless of what the value of $C_{1}$ was because
$\lim(\rho^{5/3})/\rho^{4/3}=\infty$. But this is no longer true with the
special relativistic kinetic energy. Instead, we have (this asymptotic can be
found in [19] chapter X)

$\lim_{\rho\rightarrow\infty}\frac{g\left(\alpha_{f}\rho_{f}^{1/3}\right)}{\rho^{4/3}}=6\left[\left(\frac{h}{m_{f}c}\right)^{4}\left(\frac{3}{8\pi}\right)^{4/3}\right],$
(49)

so if we have

$(m_{p}N_{p}+m_{e}N_{e})^{2/3}<\frac{\pi
2^{2/3}hc\left(\frac{3}{8\pi}\right)^{4/3}}{(GKm_{p}^{4/3})}$ (50)

Then we can apply the same procedure as in the nonrelativistic case to bound
the $||\rho_{f}||_{L^{4/3}}$ norms.

  3. 3.

Because $g(z)\geq 0$ for $z\geq 0$, we can again apply Fatou’s lemma to
conclude that

$\int g(\alpha_{f}\rho_{f}^{1/3})\leq\liminf\int
g(\alpha_{f}\rho^{n}_{f})^{1/3})$ (51)

  4. 4.

For proving that $\int\rho_{p}=N_{p}$ and $\int\rho_{e}=N_{e}$, the case that
$\int\rho_{p}<N_{p}$ and $\int\rho_{e}<N_{e}$ is exactly the same. The case
that either $\int\rho_{p}<N_{p}$ while $\int\rho_{e}=N_{e}$ or
$\int\rho_{p}=N_{p}$ and $\int\rho_{e}<N_{e}$ can be adapted using the same
arguments as given by item 1. on this list.

###  3.2 Critical mass estimate

In the special relativistic case, it is not too difficult to prove something
of a converse to Theorem 3.5. This is the existence of a mass limit analogous
to the Chandrasekhar limit found in the effective singular density case [19].

Fix $(N_{e},N_{p})$. If we take $\rho_{f}=C_{f}\chi_{B_{R}}$ for
$C_{f}=\frac{3N_{f}}{4\pi R^{3}}=:\beta_{f}^{3}/R^{3}$, then we can get an
exact expression for the energy. The electric and gravitational energies
simplify to $(3/10R)(q^{2}-GM^{2})$, where $Q$ is the total charge of the
configuration and $M$ is the total mass of the configuration. The kinetic
energies we compute as the unappealing

$\frac{4}{3}\pi R^{3}\sum_{f}\frac{\pi
m^{4}_{f}c^{5}}{3h^{3}}\left[8\alpha_{f}^{3}\frac{\beta_{f}^{3}}{R^{3}}\left(\left(\alpha_{f}^{2}\frac{\beta_{f}^{2}}{R^{2}}+1\right)^{1/2}-1\right)-\alpha_{f}\frac{\beta_{f}}{R}\left(2\alpha_{f}^{2}\frac{\beta_{f}^{2}}{R^{2}}-3\right)\left(\alpha_{f}^{2}\frac{\beta_{f}^{2}}{R^{2}}+1\right)^{1/2}-\sinh^{-1}\left(\alpha_{f}\frac{\beta_{f}}{R}\right)\right].$
(52)

It is however, easy to see that this term behaves like $K/R$ for some constant
$K$ when $R$ is very small, see again [19] chapter X for the asymptotic
expression. Therefore, if we let $R\rightarrow 0$, we see that the potential
and kinetic energies are proportional. By using that
$\lim_{z\rightarrow\infty}A(z)\sim 6z^{4}$, we can compute that

$K=\frac{3^{2/3}}{2^{10/3}}\frac{hc}{\pi^{2/3}}(N_{e}^{4/3}+N_{p}^{4/3}).$
(53)

By comparing $K$ to $(3/10)(q^{2}-GM^{2})$, we can see that, if the
gravitational potential dominates the electric and kinetic energies, letting
$R\rightarrow 0$ gives us a sequence of configurations with energy unbounded
below. Then for such combinations of $N_{e}$ and $N_{p}$, no minimizer could
exist. This result is analogous to the Chandrasekhar Mass of the single
density model:

$-\left(\frac{hc}{G}\right)^{3/2}\frac{2^{7/2}}{3^{1/2}\pi}\left(\xi^{2}\frac{d\theta_{3}}{d\xi}\right)_{\xi=\xi_{3}}\frac{1}{(\mu_{e}H)^{2}}$
(54)

where $\theta_{3}$ is the solution to the normalized Lane-Emden equation of
index 3, $\xi_{3}$ is the first zero of this solution, $\mu_{e}$ is the
average molecular weight per electron, and $H$ is the mass of the hydrogen
atom. The last fraction in particular shows how this value depends on the
composition of the star; in our case that would correspond to depending on
$N_{p}/N_{e}$. In fact, if we write $N_{p}/N_{e}=\mathfrak{r}$, we can compute
that the energy of configurations with this ratio is unbounded below if

$\left(\frac{hc}{G}\right)^{3/2}\frac{5^{3/2}}{2^{7/2}3^{1/2}\pi}\frac{(\mathfrak{r}^{4/3}+1)^{3/2}}{\left[(m_{p}\mathfrak{r}+m_{e})^{2}-\frac{C}{G}q^{2}(\mathfrak{r}-1)^{2}\right]^{3/2}}<N_{e}$
(55)

This tells us that for each admissible ratio $N_{e}/N_{p}$, like in the
Chandrasekhar case, there is some upper bound on the total mass of the
configuration, a crude estimate of which is given by the left hand side of
(55) multiplied by $(\mathfrak{r}m_{p}+m_{e})$.

##  4 Minimizers of the Energy Functional with Compact Support

This section presents a proof of the existence of radially symmetric
minimizers of (3) and its special relativistic version. It is a use of the
general technique outlined in sections 5 and 6 of [14] which prove the
existence of minimizers of functionals of the form (2), subject to certain
constraints on $A$, $\rho$, and $J$. However, there is an important difference
in the result of this section as compared with the previous. In the process of
proving existence, we also prove that there is a minimizer with compact
support. This comes at the cost of losing the equalities in (25); that is, the
technique cannot show that when the ratio $N_{e}/N_{p}$ equals either of the
bounds in (25) a compactly supported minimizer exists. This is because the
unique minimizer does not have compact support when $N_{p}/N_{e}$ saturates
either of the bounds in (25).

For $R>0$ and fixed $(N_{e},N_{p})$, let $W_{R}$ be the set of pairs
$(\rho_{e},\rho_{p})$ such that $\int\rho_{f}=N_{f}$, $\rho_{f}\geq 0$,
$\rho_{f}$ is radially symmetric, $\rho_{f}=0$ outside the ball of radius $R$,
and such that $\rho_{f}\leq R$; note that unlike in the previous section,
these functions integrate to $N_{f}$ exactly. This set will be nonempty for
large enough $R$. Define $W$ to be as $W_{R}$ without the $R$ related
constraints.

######  Lemma 4.1.

For $R$ large enough, there is a pair $(\rho^{R}_{e},\rho^{R}_{p})\in W_{R}$
which minimizes $E_{(N_{e},N_{p})}$ on $W_{R}$.

###### Proof.

It is straightforward to see that $W_{R}$ is convex, and, fixing $1<p<\infty$,
that $W_{R}$ is a closed and bounded subset of $L^{p}\times L^{p}$. Since any
bounded and closed subset of $L^{p}\times L^{p}$ is weakly compact by Banach-
Alaoglu, we show that $E_{(N_{e},N_{p})}$ is lower semicontinuous in the weak
topology for $p=\frac{5}{3}$.

Since we are considering functions on a bounded set, it is clear that if
$\rho_{e,n}$ converges weakly in $L^{5/3}$ to $\rho_{e}$, we have

$\int_{\mathbb{R}^{3}}\rho^{5/3}_{e,n}d^{3}x\rightarrow\int_{\mathbb{R}^{3}}\rho^{5/3}_{e}d^{3}x.$
(56)

The same is true for the second term in the functional.

For the third term, note that the operator

$B:\rho\rightarrow\int_{\mathbb{R}^{3}}\frac{\rho(y)}{|x-y|}d^{3}y$

is a compact operator from $L^{5/3}$ to $L^{\infty}$, by Sobolev embedding.
Suppose that $\\{(\rho_{e,n},\rho_{p,n})\\}$ converges weakly to
$(\rho_{e},\rho_{p})$. The sequence $\\{\rho_{p,n}-\rho_{e,n}\\}$ is thus
bounded. Therefore $B(\rho_{p,n}-\rho_{e,n})\rightarrow B(\rho_{p}-\rho_{e})$
in $L^{\infty}$. Then we can say

$\displaystyle\int_{\mathbb{R^{3}}}|(\rho_{p,n}-\rho_{e,n})(x)B(\rho_{p,n}-\rho_{e,n})(x)-(\rho_{p}-\rho_{e})(x)B(\rho_{p}-\rho_{e})(x)|d^{3}x$
$\displaystyle\leq\int_{\mathbb{R^{3}}}|\rho_{p,n}(x)-\rho_{e,n}(x)||B(\rho_{p,n}-\rho_{e,n})(x)-B(\rho_{p}-\rho_{e})(x)|d^{3}x$
(57)
$\displaystyle+\int_{\mathbb{R^{3}}}|(\rho_{p,n}(x)-\rho_{e,n}(x))-(\rho_{p}-\rho_{e})(x))||B(\rho_{p}-\rho_{e})(x)|d^{3}x$
$\displaystyle\rightarrow 0$

So the third and also fourth terms of the functional are even continuous in
the weak topology.

We may therefore conclude that if $\\{(\rho_{e,n},\rho_{p,n})\\}$ is a
minimizing sequence in $W_{R}$ which converges weakly to
$(\rho_{e},\rho_{p})$,

$E_{(N_{e},N_{p})}(\rho_{e},\rho_{p})\leq\liminf
E_{(N_{e},N_{p})}(\rho_{e,n},\rho_{p,n})$

so $(\rho_{e},\rho_{p})$ is a minimizer, which we now call
$(\rho^{R}_{e},\rho^{R}_{p})$. ∎

The next step, again following [14], is to show that there is some $\hat{R}$
such that for $R\geq\hat{R}$, the minimizer in $W_{R}$ lies in $W_{\hat{R}}$.
Then a limiting argument can be applied to show this function will be the
minimizer when we consider functions on the entire space. To prove this, we
use (4) and (5), but restricted to $W_{R}$, to say that

$\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)=-q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y+Gm_{p}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y+\lambda^{R}_{p}$
(58)

where $\rho_{p}^{R}>0$, and

$\frac{5}{3}k_{e}(\rho_{e}^{R})^{2/3}(x)=q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y+Gm_{e}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y+\lambda^{R}_{e}$
(59)

where $\rho_{e}^{R}>0$. We will show that, under certain conditions, the right
sides of (58) and (59) become strictly negative for large $R$, so we must have
that $\rho_{e}^{R}$ and $\rho_{p}^{R}$ vanish for large $R$.

We will also need that, on the sets where $\rho_{p}^{R}=0$ and
$\rho_{e}^{R}=0$, we have

$\lambda_{p}^{R}\leq
q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y-Gm_{p}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y$
(60)

and

$\lambda_{e}^{R}\leq-q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y-Gm_{e}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y,$
(61)

respectively. See [14] for details.

We have

######  Theorem 4.2.

If

$\frac{1-\frac{Gm_{p}^{2}}{q^{2}}}{1+\frac{Gm_{e}m_{p}}{q^{2}}}<\frac{N_{e}}{N_{p}}<\frac{1+\frac{Gm_{e}m_{p}}{q^{2}}}{1-\frac{Gm_{q}^{2}}{q^{2}}},$
(62)

there is a pair $(\rho_{e},\rho_{p})$ which minimizes (3) in $W$. Moreover,
both $\rho_{e}$ and $\rho_{p}$ have compact support.

This theorem is proved using several lemmas.

######  Lemma 4.3.

There is a constant $k_{0}$ such that
$||\rho_{e}^{R}||_{4/3},||\rho_{p}^{R}||_{4/3}\leq k_{0}$ for $R\geq R_{0}$.

###### Proof.

First note that $E(\rho_{e}^{R},\rho_{p}^{R})\leq
E(\rho_{e}^{R_{0}},\rho_{p}^{R_{0}})$ when $R\geq R_{0}$. From (3), we get two
inequalities

$k_{e}\int_{\mathbb{R}^{3}}(\rho_{e}^{R})^{5/3}(x)d^{3}x\leq\frac{G}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(x)+m_{e}\rho^{R}_{e}(x))(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}yd^{3}x+E(\rho_{e}^{R_{0}},\rho_{p}^{R_{0}})$
(63)

and

$k_{p}\int_{\mathbb{R}^{3}}(\rho_{p}^{R})^{5/3}(x)d^{3}x\leq\frac{G}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(x)+m_{e}\rho^{R}_{e}(x))(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}yd^{3}x+E(\rho_{e}^{R_{0}},\rho_{p}^{R_{0}})$
(64)

We apply proposition 6 of [14] to conclude

$\displaystyle\frac{G}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(x)+m_{e}\rho^{R}_{e}(x))(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}yd^{3}x$
$\displaystyle\leq\hat{c}\left(\int_{\mathbb{R^{3}}}|m_{p}\rho^{R}_{p}(x)+m_{e}\rho^{R}_{e}(x)|^{4/3}d^{3}x\right)(m_{e}N_{e}+m_{p}N_{p})^{2/3}$
(65)
$\displaystyle\leq\hat{c}||\rho_{e}^{R}||^{4/3}_{4/3}+\hat{c}||\rho_{p}^{R}||^{4/3}_{4/3}$
(66)

where we have absorbed constants into $\hat{c}$. There are three cases to
consider. In the first, both $||\rho_{e}^{R}||_{4/3}$ and
$||\rho_{p}^{R}||_{4/3}$ are bounded as $R\rightarrow\infty$, so there is
nothing to prove. In the second, one of $||\rho_{e}^{R}||_{4/3}$ or
$||\rho_{p}^{R}||_{4/3}$ is bounded while the other is unbounded. In the third
case, both are unbounded.

Let us consider the second case, and suppose $||\rho_{e}^{R}||_{4/3}$ is
unbounded. For any $\hat{c}$, there is some value $s_{\epsilon}$ such that

$k_{e}s^{5/3}\geq(\hat{c}+2\epsilon)s^{4/3}$

when $s>s_{\epsilon}$. Therefore,

$\displaystyle\int_{\mathbb{R}^{3}}|\rho_{e}^{R}(x)|^{4/3}d^{3}x$
$\displaystyle=\int_{\rho_{e}^{R}<s_{\epsilon}}|\rho_{e}^{R}(x)|^{4/3}d^{3}x+\int_{\rho_{e}^{R}>s_{\epsilon}}|\rho_{e}^{R}(x)|^{4/3}d^{3}x$
(67) $\displaystyle\leq
N_{e}s_{\epsilon}^{1/3}+(\hat{c}+2\epsilon)^{-1}k_{e}\int_{\mathbb{R}^{3}}|\rho_{e}^{R}(x)|^{5/3}d^{3}x$
(68)

For $||\rho_{e}^{R}||_{4/3}^{4/3}$ sufficiently large, this implies (note that
$2\epsilon$ is replaced with $\epsilon$)

$(\hat{c}+\epsilon)\int_{\mathbb{R}^{3}}|\rho_{e}^{R}(x)|^{4/3}d^{3}x\leq
k_{e}\int_{\mathbb{R}^{3}}|\rho_{e}^{R}(x)|^{5/3}d^{3}x$ (69)

Plugging (69) into (63), we get

$\epsilon||\rho_{e}^{R}||_{4/3}^{4/3}\leq
E(\rho_{e}^{R_{0}},\rho_{p}^{R_{0}})+\hat{c}||\rho_{p}^{R}||_{4/3}^{4/3}$ (70)

which is a contradiction if $||\rho_{e}^{R}||_{4/3}$ is unbounded and
$||\rho_{p}^{R}||_{4/3}$ is bounded. The case with the roles of the densities
reversed is of course the same.

Now let us consider the case when both $||\rho_{p}^{R}||_{4/3}$ and
$||\rho_{e}^{R}||_{4/3}$ are unbounded. We will apply the same sort of ideas
as in the previous case, but slightly modified. Since both norms are
unbounded, we can use the same argument to show that for $R$ large enough,

$\left(\hat{c}+\frac{\hat{c}}{1-\epsilon}\right)\int_{\mathbb{R}^{3}}|\rho_{e}^{R}(x)|^{4/3}d^{3}x\leq
k_{e}\int_{\mathbb{R}^{3}}|\rho_{e}^{R}(x)|^{5/3}d^{3}x$ (71)

and

$2\hat{c}\int_{\mathbb{R}^{3}}|\rho_{p}^{R}(x)|^{4/3}d^{3}x\leq
k_{p}\int_{\mathbb{R}^{3}}|\rho_{p}^{R}(x)|^{5/3}d^{3}x$ (72)

Using these inequalities with (63) and (64), we get

$\frac{\hat{c}}{1-\epsilon}||\rho_{e}^{R}||_{4/3}^{4/3}\leq
E(\rho_{e}^{R_{0}},\rho_{p}^{R_{0}})+\hat{c}||\rho_{p}^{R}||_{4/3}^{4/3}$ (73)

and

$\hat{c}||\rho_{p}^{R}||_{4/3}^{4/3}\leq
E(\rho_{e}^{R_{0}},\rho_{p}^{R_{0}})+\hat{c}||\rho_{e}^{R}||_{4/3}^{4/3}$ (74)

which again gives a contradiction. This completes the proof. ∎

######  Lemma 4.4.

$\lambda_{e}^{R}<0$ and $\lambda_{p}^{R}<0$ for all $R$.

###### Proof.

Fix $R$ and consider $|x|=2R$. Since $\rho_{e}^{R}(x)=\rho_{p}^{R}(x)=0$, we
have that $\lambda_{p}^{R}$ and $\lambda_{e}^{R}$ satisfy (60) and (61),
respectively. Since our densities have spherical symmetry, we can apply
Newton’s theorem to these integrals to get

$\lambda_{p}^{R}\leq
q^{2}\frac{N_{p}-N_{e}}{2R}-Gm_{p}\frac{m_{p}N_{p}+m_{e}N_{e}}{2R}$ (75)

and

$\lambda_{e}^{R}\leq-q^{2}\frac{N_{p}-N_{e}}{2R}-Gm_{e}\frac{m_{p}N_{p}+m_{e}N_{e}}{2R}$
(76)

The result follows from (25). ∎

######  Lemma 4.5.

$||\rho_{e}^{R}||_{\infty},||\rho_{p}^{R}||_{\infty}<k_{1}$ for all $R>R_{0}$.

###### Proof.

We combine (58) and (59) with the results of Lemma (4.4). Since
$\lambda_{e}^{R}<0$ and $\lambda_{p}^{R}<0$, we have

$\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)\leq\max\\{0,-q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y+Gm_{p}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y\\}$
(77)

and

$\frac{5}{3}k_{e}(\rho_{e}^{R})^{2/3}(x)\leq\max\\{0,q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y+Gm_{e}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y\\}$
(78)

as the densities cannot be negative. Since we have shown that
$||\rho_{f}||_{4/3}$ is bounded, we apply proposition 5 with $p=4/3$ from [14]
to the right hand side of (77) to conclude (choosing $r=6$ for concreteness)

$\displaystyle||B(-q^{2}(\rho_{p}^{R}-\rho_{p}^{R})+Gm_{p}(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R}))||_{6}$
$\displaystyle\leq
c_{0}||(-q^{2}(\rho_{p}^{R}-\rho_{p}^{R})+Gm_{p}(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})||_{1}^{b}||(-q^{2}(\rho_{p}^{R}-\rho_{p}^{R})+Gm_{p}(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})||_{4/3}^{(1-b)}$
(79)
$\displaystyle+c_{0}||(-q^{2}(\rho_{p}^{R}-\rho_{p}^{R})+Gm_{p}(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})||_{1}^{c}||(-q^{2}(\rho_{p}^{R}-\rho_{p}^{R})+Gm_{p}(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})||_{4/3}^{1-c}$
$\displaystyle\leq
c_{0}\left[q^{2}(N_{p}+N_{e})+Gm_{p}(m_{p}N_{P}+m_{e}N_{e})\right]^{b}\left[2q^{2}k_{0}^{2}+2Gm_{p}^{2}k_{0}^{2}\right]^{1-b}$
(80)
$\displaystyle+c_{0}\left[q^{2}(N_{p}+N_{e})+Gm_{p}(m_{p}N_{P}+m_{e}N_{e})\right]^{c}\left[2q^{2}k_{0}^{2}+2Gm_{p}^{2}k_{0}^{2}\right]^{1-c}.$

So we can conclude that the $L^{6}$ norms of these potentials are uniformly
bounded for all $R$ large enough; the same result holds for (78).

From here we use a bootstrapping procedure. The results of the previous
paragraph combined with (77) and (78) imply that both $\rho_{p}^{R}$ and
$\rho_{e}^{R}$ are uniformly bounded in the $L^{4}$ norm. We then apply the
second conclusion of proposition 5 in [14] and the same inequalities as above
but with $4/3$ replaced with $4$ to conclude that the right hand sides of (77)
and (78) are uniformly bounded for all $R$ in $L^{\infty}$ norm. But this
implies that $||\rho_{e}^{R}||_{\infty}$ and $||\rho_{p}^{R}||_{\infty}$ are
also both uniformly bounded. ∎

######  Lemma 4.6.

There is an $e<0$ such that $E(\rho_{e}^{R},\rho_{p}^{R})\leq e$ for all large
$R$.

###### Proof.

It is easy to see that this is true in the case that $N_{e}=N_{p}$. For if we
take any configuration in $W$ such that $\rho_{e}=\rho_{p}$, the electric
energy is zero and we are left with

$E(\rho_{e},\rho_{p})=(k_{e}+k_{p})\int_{\mathbb{R}^{3}}\rho^{5/3}_{e}(x)d^{3}x-\frac{G}{2}(m_{e}+m_{p})^{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{\rho_{e}(x)\rho_{e}(y)}{|x-y|}d^{3}yd^{3}x$
(81)

Using the same “spreading out” argument we used when we considered the zero-
gravity case, we see that as $\lambda\rightarrow\infty$, the kinetic energy
terms go to zero like $1/\lambda^{2}$ while the gravitational energy goes to
zero like $1/\lambda$. Then for large enough $\lambda$,
$E(\rho_{e}^{\lambda},\rho_{p}^{\lambda})<0$. Since this procedure works for
any neutral configuration, it works for one of compact support, and so we must
take $R$ larger than this $\lambda$ if we want to find a configuration with
negative energy in $W_{R}$.

Therefore, because each of the terms in (3) is in some sense continuous in
$(\rho_{e},\rho_{p})$ (without trying to make this precise in any way, we just
use the intuition), there is some interval of values for $N_{e}/N_{p}$
containing 1 such that $\mathcal{E}<0$. We determine the size of this interval
by finding example functions in $W$.

Suppose that for a fixed $(N_{e},N_{p})$ such that the ratio $N_{e}/N_{p}$
satisfies $(\ref{5/3bounds})$, we have a pair $(\rho_{e},\rho_{p})$ both with
compact support such that $E(\rho_{e},\rho_{p})<0$. We claim we can add a
small positive or negative charge outside of this configuration and the energy
will not increase. Formally, the idea is to see how many “test particles” we
can bring in from infinity and still have negative energy. Of course, this is
not a rigorous argument, so we now make this idea rigorous.

Let us assume we are adding positive charge. Since our space $W_{R}$ consists
of spherically symmetric functions, we must add our charge in a spherically
symmetric configuration. The configuration to which we are adding has compact
support, suppose it is contained in a ball of radius $R_{0}$. Then add to
$(\rho_{e},\rho_{p})$ a particle density described by the function
$h=\epsilon\chi_{S}$ for $S=\\{x|R_{0}<|x|<R_{0}+\eta\\}$. Then we compute

$\displaystyle E(\rho_{e},\rho_{p}+g)-E(\rho_{e},\rho_{p})$
$\displaystyle=k_{p}\epsilon^{5/3}\frac{4}{3}\pi\left[(R_{0}+\eta)^{3}-R_{0}^{3}\right]+\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}\int_{S}\int_{S}\frac{1}{|x-y|}d^{3}yd^{3}x$
$\displaystyle+\epsilon\int_{S}\int_{\mathbb{R}^{3}}\frac{q^{2}(\rho_{p}-\rho_{e})-Gm_{p}(m_{p}\rho_{p}+m_{e}\rho_{e})}{|x-y|}$
(82)
$\displaystyle=k_{p}\epsilon^{5/3}\frac{4}{3}\pi\left[3R_{0}^{2}\eta+3R_{0}\eta^{2}+\eta^{3}\right]$
$\displaystyle+\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}4\pi\int_{R_{0}}^{R_{0}+\eta}r^{2}\left[\frac{4}{3}\frac{\pi}{r}\left[r^{3}-R_{0}^{3}\right]+2\pi\left[(R_{0}+\eta)^{2}-r^{2}\right]\right]dr$
$\displaystyle+4\pi\epsilon\int_{R_{0}}^{R_{0}+\eta}\left[q^{2}(N_{p}-N_{e})-Gm_{p}(m_{p}N_{p}+m_{e}N_{e})\right]rdr$
(83)
$\displaystyle=k_{p}\epsilon^{5/3}\frac{4}{3}\pi\left[3R_{0}^{2}\eta+3R_{0}\eta^{2}+\eta^{3}\right]$
$\displaystyle+\left[\frac{q^{2}}{2}-\frac{G}{2}\right]\epsilon^{2}4\pi\left[4\pi
R_{0}^{3}\eta^{2}+\frac{16\pi}{3}R_{0}^{2}\eta^{3}+\frac{8\pi}{3}R_{0}\eta^{4}+\frac{3\pi}{5}\eta^{5}\right]$
$\displaystyle-4\pi\epsilon
N_{e}(q^{2}-Gm_{p}^{2})\left(B_{p}-\frac{N_{p}}{N_{e}}\right)(2R_{0}\eta+\eta^{2})$
(84)

where $B_{p}$ is the multiplicative inverse of the left side of (25). If we
assume $\left(B_{p}-\frac{N_{p}}{N_{e}}\right)>0$, the last term is strictly
negative. Therefore, we may take $\epsilon>0$ small enough so that the
difference is negative. So we can add some positive amount of charge to the
configuration and the energy will become more negative. The same will be true
of adding a small negative amount of charge using the strict version of the
right half of (25).

We intend to treat this as an iterative procedure, and note that if
$N_{e}=N_{p}$, we can assume that the densities have compact support; the size
of the support makes no difference for the argument that
$\mathcal{E}^{N_{e},N_{e}}<0$. The issue remaining is whether we can saturate
(25) using this procedure.

To show that we can, we let $\eta=\epsilon$. Then if we have only added a
finite amount of charge to the configuration, we can also be assured the
resulting configuration has finite radius. So no matter how much charge we
add, we can always apply the procedure again as long as the strict version of
(25) is satisfied. This will give us (25). So we need to show that the
procedure can be carried out if $\eta=\epsilon$, which is to say that we can
find an $\epsilon$ so that in this case (4) is negative.

The leading order in $\epsilon$ parts of the three terms of (4) behave as
$\epsilon^{8/3}R_{0}^{2}$, $\epsilon^{3}R_{0}^{4}$, and $\epsilon^{2}R_{0}$,
respectively. So if we take $\epsilon$ to be $o(R_{0}^{-4})$, we can make the
first two terms smaller in absolute value than the last, negative term. ∎

######  Lemma 4.7.

$||B(m_{e}\rho_{e}^{R}+m_{p}\rho_{p}^{R})||_{\infty}\geq-\frac{2e}{G(m_{p}N_{p}+m_{e}N_{p})}$
for all large $R$.

###### Proof.

From the previous lemma, we have

$\displaystyle-e$
$\displaystyle\leq-E(\rho_{e}^{R},\rho_{p}^{R})\leq\frac{G}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(x)+m_{e}\rho^{R}_{e}(x))(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}yd^{3}x$
(85)
$\displaystyle\leq\frac{G}{2}(m_{p}N_{p}+m_{e}N_{p})||B(m_{e}\rho_{e}^{R}+m_{p}\rho_{p}^{R})||_{\infty}$
(86)

where in the second line we have applied Holder’s inequality. ∎

######  Lemma 4.8.

There is an $\epsilon>0$ such that for large $R$,

$\int_{|x-x_{R}|<1}\rho_{f}^{R}(x)d^{3}x\geq\epsilon$ (87)

for $f=e$ or $p$ and some $x_{R}$.

###### Proof.

The point of this lemma is to show that the densities cannot spread out
indefinitely; at least one must have some mass accumulated somewhere.

Define

$\epsilon_{e}^{R}=\sup_{x}\int_{|y-x|<1}\rho_{e}^{R}(y)d^{3}y$ (88)

and $\epsilon_{p}^{R}$ analogously. We have

$\displaystyle B(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(x)$
$\displaystyle=\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y$
(89)
$\displaystyle=\int_{|x-y|<1}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y+\int_{1<|x-y|<r}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y$
(90)
$\displaystyle+\int_{|x-y|>r}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y$

The last term is bounded by $\frac{m_{p}N_{p}+m_{e}N_{e}}{r}$. The set
$1<|x-y|<r$ can be covered by $cr^{3}$ balls of radius 1, and can therefore be
bounded by $cr^{3}(m_{p}\epsilon_{p}^{R}+m_{e}\epsilon_{e}^{R})$. Since lemma
(4.5) shows that $\rho_{e}^{R}$ and $\rho_{p}^{R}$ are bounded in the
$L^{\infty}$ norm, we can apply proposition 5 with $p=r=\infty$ of [14] on the
ball of radius 1 to get

$\int_{|x-y|<1}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y\leq
k_{4}((m_{e}\epsilon_{e}^{R}+m_{p}\epsilon_{p}^{R})^{b}+(m_{e}\epsilon_{e}^{R}+m_{p}\epsilon_{p}^{R})^{c})$
(91)

Combining these three estimates, we get

$B(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(x)\leq\frac{m_{p}N_{p}+m_{e}N_{e}}{r}+cr^{3}(m_{p}\epsilon_{p}^{r}+m_{e}\epsilon_{e}^{R})+k_{4}((m_{e}\epsilon_{e}^{R}+m_{p}\epsilon_{p}^{R})^{b}+(m_{e}\epsilon_{e}^{R}+m_{p}\epsilon_{p}^{R})^{c})$
(92)

Now if both $\epsilon_{e}^{R}$ and $\epsilon_{p}^{R}$ could get arbitrarily
close to 0 as $R$ became large, then we could choose $r$ to be very large, and
then $\epsilon_{e}^{R}$ and $\epsilon_{p}^{R}$ to be very small and we would
get that

$B(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(x)<-\frac{2e}{G(m_{p}N_{p}+m_{e}N_{p})}$
(93)

contradicting the previous lemma.

Note that because of spherical symmetry, we also get a bound on $|x_{R}|$ for
all $R$; we will say that $|x_{R}|<r_{0}$. ∎

######  Lemma 4.9.

For some $l<0$, $\lambda_{e}^{R}+\lambda_{p}^{R}<l$ for large $R$.

###### Proof.

Suppose $R>>r_{0}$. If $|x-x_{R}|<r<R$, then we can say that

$\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y\geq\int_{|y-x_{R}|<1}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y\geq\frac{m_{e}\epsilon}{r+1}$
(94)

which is independent of $R$, but of course depends on $x$. As $|x-x_{R}|$
increases (and $r$ with it), we see that this value decreases as $r^{-1}$.
Because $\int\rho^{R}_{p}=N_{p}$ and $\int\rho^{R}_{e}=N_{e}$, in the
$r$-sphere around $x_{R}$ there must be some point $x$ such that
$\rho_{f}^{R}(x)<kr^{-3}$, for $k$ independent of $R$. So if we take $r$ large
enough, we can find an $x$ such that the sum of the kinetic energy terms is
less than half the negative of the gravitational term.

We want to find the same $x$ for both $\rho^{R}_{e}$ and $\rho^{R}_{p}$, so we
prove that for any $r$ large enough, there is always a set $U_{f}$ such that
$m(U_{f})>\frac{2}{3}\pi r^{3}$, where $m$ means the measure, and on $U_{f}$,
$\rho_{f}^{R}<k_{1}r^{-3}$. Then $m(U_{p}\cap U_{e})\neq 0$, so there must be
some $\overline{x}$ such that $\rho_{e}^{R}(\overline{x})<k_{1}r^{-3}$ and
$\rho_{p}^{R}(\overline{x})<k_{1}r^{-3}$. Assume this is not true. Then there
is a set $S_{f}$ such that $m(S_{f})\geq\frac{2}{3}\pi r^{3}$ and
$\rho_{f}^{R}>k_{1}r^{-3}$ on $S_{f}$. For the correctly chosen $k_{1}$ (still
independent of $R$) this is a contradiction.

In the “worst case” scenario, at $\overline{x}$ both densities are positive
and we add equations (58) and (59) together to obtain

$\lambda_{e}^{R}+\lambda_{p}^{R}=\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(\overline{x})+\frac{5}{3}k_{e}(\rho_{e}^{R})^{2/3}(\overline{x})-G(m_{e}+m_{p})\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|\overline{x}-y|}d^{3}y$
(95)

If one of the densities is 0, then (95) is not true. But this would mean that
$|x|>|\overline{x}|$ implies $\rho_{f}^{R}(x)=0$. As we intend to use (95) to
show that when $R$ is large one of the densities must be zero for $|x|$ large,
this is no problem.

So for any $R$ large enough, we have found an $\overline{x}$ such that the
right side is negative, and therefore $\lambda_{e}^{R}+\lambda_{p}^{R}<l<0$. ∎

######  Lemma 4.10.

For $R$ large enough, there is an $\alpha$ such that for $|x|>\alpha$,
$B(m_{p}\rho^{R}_{p}+m_{e}\rho^{R}_{e})(x)\leq-l$.

###### Proof.

Define $\epsilon_{\alpha}^{R}$ for $\alpha<<R$ by

$\epsilon_{\alpha}^{R}=\sup_{|x|>\alpha/2}\int_{|x-y|<1}(m_{e}\rho_{e}^{R}+m_{p}\rho_{p}^{R})(y)dy$
(96)

Radial symmetry implies that $\epsilon^{R}_{\alpha}$ decays faster than
$1/\alpha$ as $\alpha$ gets larger. Then if we take $|x|>\alpha$, using the
same arguments as in lemma 4.8, we find

$B(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(x)\leq
c_{1}((\epsilon_{\alpha}^{R})^{a}+(\epsilon_{\alpha}^{R})^{b})+c_{1}\epsilon_{\alpha}^{R}r^{3}+c_{3}r^{-1}$
(97)

So we can take first $r$ very large, and then $\alpha$ very large to get that
$B(m_{p}\rho^{R}_{p}+m_{e}\rho^{R}_{e})(x)\leq-l$, which holds for all $R$
large enough. ∎

As long as both densities are nonzero, we can add (58) and (59) together and
rewrite it as

$\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)+\frac{5}{3}k_{e}(\rho_{e}^{R})^{2/3}(x)=G(m_{e}+m_{p})\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho^{R}_{p}(y)+m_{e}\rho^{R}_{e}(y))}{|x-y|}d^{3}y+\lambda_{e}^{R}+\lambda_{p}^{R}$
(98)

But by the last lemma, when $R$ is very large, for $|x|>\alpha$ the left side
would be negative, so one of the densities must become zero at a smaller $|x|$
so that this is no longer an equality. Then for $|x|$ large enough, we may
assume that one of the densities is zero, although we do not know which one.
The final step is to prove that there is a $\beta$ such that if $|x|>\beta$,
both densities are 0.

We need

######  Lemma 4.11.

Both densities have contiguous support.

###### Proof.

Suppose that $\rho_{e}^{R}$ did not have contiguous support. By this we mean
that there is some open set $A$ such that $\rho_{e}^{R}=0$ on $A$, but that
there is another open set $O$ on which $\rho_{e}^{R}>0$ and such that

$\sup\\{|x|:x\in A\\}\leq\inf\\{|x|:x\in O\\}$ (99)

Suppose $O$ is such that this is actually an equality ($O$ and $A$ are
adjacent).

The idea is that we can move some amount of $\rho_{e}^{R}$ and lower the total
energy contradicting that $(\rho_{e}^{R},\rho_{p}^{R})$ is a minimizer.
Consider

$v_{1}(x)=\begin{cases}-\rho_{e}^{R}&x\in K_{O}\\\ \rho_{e}^{R}(x+\alpha)&x\in
K_{A}\\\ 0&\text{else}\end{cases}$

where $K_{A}$ is some small set contained in $A$, and $K_{O}$ is $K_{A}$
shifted by the amount $\alpha$ so that it lies in $O$. Then consider

$\displaystyle\delta_{1}$
$\displaystyle=E(\rho_{e}^{R}+v_{1},\rho_{p}^{R})-E(\rho_{e}^{R},\rho_{p}^{R})$
(100)
$\displaystyle=-q^{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{1}(x)(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))}{|x-y|}d^{3}y-Gm_{e}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{1}(x)(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y))}{|x-y|}d^{3}y$
(101)
$\displaystyle+\frac{q^{2}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{1}(x)v_{1}(y)}{|x-y|}d^{3}y-\frac{Gm^{2}_{e}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{1}(x)v_{1}(y)}{|x-y|}d^{3}y$
$\displaystyle=\int_{K_{O}}\rho_{e}^{R}(x)\int_{\mathbb{R}^{3}}\frac{q^{2}(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))+Gm_{e}(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y))}{|x-y|}d^{3}y$
(102)
$\displaystyle-\int_{K_{A}}\rho_{e}^{R}(x+\alpha)\int_{\mathbb{R}^{3}}\frac{q^{2}(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))+Gm_{e}(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y))}{|x-y|}d^{3}y$
$\displaystyle+\frac{q^{2}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{1}(x)v_{1}(y)}{|x-y|}d^{3}y-\frac{Gm^{2}_{e}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{1}(x)v_{1}(y)}{|x-y|}d^{3}y$

As we can make the sets $K_{A}$ and $K_{O}$ as small as we like, the first two
terms in the last equality can be made to dominate the last two, so let us
focus on the first two. Now if $\delta_{1}<0$ for any choice of $K_{A}$ and
$K_{O}$, then $(\rho_{e}^{R},\rho_{p}^{R})$ is not a minimizer after all. So
assume that $\delta_{1}\geq 0$.

We now split the proof into two cases. In the first, there is not an open set
$D$ on which $\rho_{e}^{R}>0$ and such that

$\sup\\{|x|:x\in D\\}\leq\inf\\{|x|:x\in A\\}$ (104)

This would imply that the $\rho_{e}^{R}=0$ at the origin. In case 2, there is
such a $D$. We handle the cases separately.

Define
$V(x)=-B(q^{2}(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))+Gm_{e}(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y)))$,
so that $V$ is the potential felt by the electron density. In case 1,
$\delta_{1}>0$ implies that $V(x)$ is smaller on $O$ than on $A$. But

$-\Delta V=-q^{2}\rho_{p}^{R}(y)-Gm_{e}m_{p}\rho_{p}^{R}(y)\leq 0$ (105)

implies that $\Delta V\geq 0$ on $A$, or that the maximum of $V$ on $A$ occurs
at its boundary. But $O$ is the boundary of $A$, so we get a contradiction.

In case 2, consider

$v_{2}(x)=\begin{cases}-\rho_{e}^{R}&x\in K_{D}\\\ \rho_{e}^{R}(x-\beta)&x\in
K_{A}\\\ 0&\text{else}\end{cases}$

where $K_{D}$ is $K_{A}$ shifted by $-\beta$ and on which $\rho_{e}^{R}>0$.
Further assume $D$ is on the boundary of $A$ (so that $O$ and $D$ are on the
boundaries of $A$). Then consider

$\displaystyle\delta_{2}$
$\displaystyle=E(\rho_{e}^{R}+v_{2},\rho_{p}^{R})-E(\rho_{e}^{R},\rho_{p}^{R})$
(106)
$\displaystyle=-q^{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{2}(x)(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))}{|x-y|}d^{3}y-Gm_{e}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{2}(x)(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y))}{|x-y|}d^{3}y$
(107)
$\displaystyle+\frac{q^{2}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{2}(x)v_{2}(y)}{|x-y|}d^{3}y-\frac{Gm^{2}_{e}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{2}(x)v_{2}(y)}{|x-y|}d^{3}y$
$\displaystyle=\int_{K_{D}}\rho_{e}^{R}(x)\int_{\mathbb{R}^{3}}\frac{q^{2}(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))+Gm_{e}(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y))}{|x-y|}d^{3}y$
(108)
$\displaystyle-\int_{K_{A}}\rho_{e}^{R}(x-\beta)\int_{\mathbb{R}^{3}}\frac{q^{2}(\rho_{p}^{R}(y)-\rho_{e}^{R}(y))+Gm_{e}(m_{p}\rho_{p}^{R}(y)+m_{e}\rho_{e}^{R}(y))}{|x-y|}d^{3}y$
$\displaystyle+\frac{q^{2}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{2}(x)v_{2}(y)}{|x-y|}d^{3}y-\frac{Gm^{2}_{e}}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{v_{2}(x)v_{2}(y)}{|x-y|}d^{3}y$

Again, we have that if $\delta_{2}<0$, then $(\rho_{e}^{R},\rho_{p}^{R})$ is
not a minimizer. So assume that $\delta_{2}\geq 0$. If $\delta_{1}>0$ and
$\delta_{2}>0$, then $V(x)$ is lesser on both $K_{O}$ and $K_{D}$ than on
$K_{A}$. But again we have that $\Delta V\geq 0$ in $A$, or that the maximum
for $V$ in $A$ is on the boundary. This is again a contradiction and one of
$\delta_{1}$ or $\delta_{2}$ must be negative. Therefore the densities must be
contiguous.

In the case that $\delta_{1}$ or $\delta_{2}$ are equal to zero,
$(\rho_{e}^{R}+v,\rho_{p}^{R})$ is also a minimizer, but one where we have
moved some electron density into $Z_{n}$. If $\rho_{e}^{R}$ is still not
contiguous, we apply the procedure again and try to find where $\delta_{1}$
and $\delta_{2}$ are both greater than zero. So we repeat until either we find
such $\delta_{1}$ and $\delta_{2}$, or have a contiguous density.

The case in which $\rho_{p}^{R}=0$ is similar. ∎

So not only can we assume that for $|x|$ large enough one of the densities is
zero, we must have that it is the same density as we increase $|x|$.

Suppose it is $\rho_{e}^{R}$, and let $L_{R}=\sup\\{|x||\rho_{e}^{R}(x)>0\\}$.
Now if we take $|x|>L_{R}$, then we find that (assuming $\rho_{p}^{R}\neq 0$)

$\displaystyle\lambda_{p}^{R}$
$\displaystyle=\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)+q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}^{R}-\rho_{e}^{R})(y)}{|x-y|}d^{3}y-Gm_{p}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}^{R}+m_{e}\rho_{e}^{R})(y)}{|x-y|}d^{3}y$
(110)
$\displaystyle=\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)+\frac{q^{2}(N_{p}(|x|)-N_{e})}{|x|}-\frac{Gm_{p}(m_{p}N_{p}(|x|)+m_{e}N_{e})}{|x|}$
(111) $\displaystyle+4\pi[q^{2}-Gm_{p}^{2}]\int_{|x|}^{R}\rho^{R}_{p}(r)rdr$
$\displaystyle\leq\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)+\frac{q^{2}(N_{p}-N_{e})}{|x|}-\frac{Gm_{p}(m_{p}N_{p}+m_{e}N_{e})}{|x|}$
(112)

where by $N_{p}(|x|)$ we mean the number of positively charged particles
within a radius $R=|x|$, and we have abused notation to define
$\rho_{p}^{R}(r)$ in the obvious way.

Since we have fixed $N_{p}$ and $N_{e}$ to satisfy the strict inequality
version of (25), we can conclude that there is some $\gamma$ such that

$q^{2}(N_{p}-N_{e})-Gm_{p}(m_{p}N_{p}+M_{e}N_{e})<\gamma<0$ (113)

Then we have

$\lambda_{p}^{R}<\frac{5}{3}k_{p}(\rho_{p}^{R})^{2/3}(x)+\frac{\gamma}{|x|}.$
(114)

Again, we can find $\overline{x}$ so that
$\rho_{p}^{R}(\overline{x})<k_{1}R^{-3}$. Therefore, we can take $R$ large
enough so we can find $\overline{x}$ showing that $\lambda_{p}^{R}<\hat{l}<0$
for some $\hat{l}$ and large $R$.

Then where $\rho_{p}^{R}>0$ we have

$\frac{q^{2}(N_{p}(|x|)-N_{e})}{|x|}-\frac{Gm_{p}(m_{p}N_{p}(|x|)+m_{e}N_{e})}{|x|}<\lambda_{p}^{R}<\hat{l}<0$
(115)

Clearly, this cannot be satisfied for large $|x|$. So for very large $R$, we
must have both densities equal to zero. The same could be concluded in the
case that $\rho_{p}^{R}$ becomes zero first.

Now the standard argument found in [14] shows that for any admissible pair
$(N_{e},N_{p})$, their exists at least one minimizer of $E$ in $W$ with
compact support.

###  4.1 The Special Relativistic Case

Just as in the first proof of the existence of minimizers of $E$, this proof
can be easily modified to accommodate the relativistic kinetic energy (12).
Again, the proof only directly uses the kinetic energy in a few places: we
address those here.

  1. 1.

In lemma 4.1, we used that $x^{5/3}$ is convex in proving lower
semicontinuity. $g(z)$ is also convex, so this step follows exactly as before.

  2. 2.

In lemma 4.3, we used that for any $\hat{c}$,
$k_{e}s^{5/3}\geq(\hat{c}+2\epsilon)s^{4/3}$ for $s>s_{\epsilon}$. Just as in
the modifications to the first proof, this must be changed. In the proof,
$\hat{c}=CG2^{1/3}m_{p}^{4/3}(m_{e}N_{e}+m_{p}N_{p})^{2/3}$, were $C$ is the
best constant in Proposition 6 of [14]. For this step of the proof to go
through, we need

$\hat{c}<\frac{\pi
m_{f}c^{5}}{3h^{3}}\lim_{s\rightarrow\infty}\frac{(g(\alpha_{f}s^{1/3}))}{s^{4/3}}$

which is true if

$(m_{e}N_{e}+m_{p}N_{p})^{2/3}<\frac{\pi
2^{2/3}hc\left(\frac{3}{8\pi}\right)^{4/3}}{GKm_{p}^{4/3}}$ (116)

This is the exact same bound as in the statement of theorem 3.5.

  3. 3.

There are a few times in the remainder of the proof where we bound
$\frac{d}{d\rho}(\rho^{5/3})$ by bounding $\rho$ by small numbers. These steps
are also the same as for small $z$, $g(z)$ behaves like $z^{5/3}$, so similar
estimates hold.

Then we have the theorem

######  Theorem 4.12.

If

$\frac{1-\frac{Gm_{p}^{2}}{q^{2}}}{1+\frac{Gm_{e}m_{p}}{q^{2}}}<\frac{N_{e}}{N_{p}}<\frac{1+\frac{Gm_{e}m_{p}}{q^{2}}}{1-\frac{Gm_{q}^{2}}{q^{2}}},$
(117)

there is a pair $(\rho_{e},\rho_{p})$ which minimizes $E^{s}$ in $W$.
Moreover, both $\rho_{e}$ and $\rho_{p}$ have compact support.

##  5 Minimizers in the Nonrelativistic Case

Let us now turn to the study of the equations (8) and (9). Based on the
analysis of the “6/3” model presented in [9], we expect there will be two
problems we will need to study. The first, following the terminology
introduced in [9], we call the bulk problem. For this problem, both densities
are positive, and we have the system given by (8) and (9). The second problem,
which we call the atmospheric problem, occurs when one of the densities is 0.
Depending on which density is zero, we ignore one of (8) or (9), and use
$\rho_{f}=0$ appropriately with the equation we focus on.

###  5.1 The Bulk Problem

With spherical symmetry and the substitution $u_{f}^{2/3}=\rho_{f}$, we can
rewrite (8) and (9) as

$u_{p}^{\prime\prime}(r)=-\frac{2}{r}u_{p}^{\prime}+\frac{12\pi}{5k_{p}}\left(q^{2}-Gm_{p}^{2}\right)u_{p}^{3/2}(r)-\frac{12\pi}{5k_{p}}\left(q^{2}+Gm_{p}m_{e}\right)u_{e}^{3/2}(r)$
(118)

and

$u_{e}^{\prime\prime}(r)=-\frac{2}{r}u_{e}^{\prime}-\frac{12\pi}{5k_{p}}\left(q^{2}+Gm_{p}m_{e}\right)u_{p}^{3/2}(r)+\frac{12\pi}{5k_{p}}\left(q^{2}-Gm^{2}_{e}\right)u_{e}^{3/2}(r).$
(119)

We want to solve this system with initial conditions

$u_{p}(0)=\alpha>0,u_{e}(0)=\beta>0,\text{ and
}u^{\prime}_{e}(0)=u_{p}^{\prime}(0)=0;$ (120)

the conditions on the derivatives necessary so we do not get a singularity at
the origin. We will first prove the well-posedness of this initial value
problem.

######  Lemma 5.1.

While both $u_{e}$ and $u_{p}$ are greater than zero, the system given by
(118), (119), and (120) is locally well-posed.

###### Proof.

Here we follow the approach of [32]. Fix $(\alpha,\beta)$. We will show that
for any positive $\alpha$ and $\beta$, there is a $\delta>0$ such that $u_{e}$
and $u_{p}$ are at least $C^{2}$ on $[0,\delta]$, and depend continuously on
the initial data. Let

$X=\\{(f,g)\in C([0,\delta])\times
C([0,\delta])|\max(||f-\alpha||_{\infty},||g-\beta||_{\infty})<\epsilon\\},$

where $\epsilon$ is chosen small enough so that $f,g>0$. Then consider
$T:X\rightarrow X$ defined as

$\displaystyle T(f,g)$
$\displaystyle=(\hat{f}(r),\hat{g}(r))=\left(\alpha-\int_{0}^{r}\int_{0}^{t}\frac{Ff^{3/2}(s)-Eg^{3/2}(s)}{t^{2}}s^{2}dsdt,\beta-\int_{0}^{r}\int_{0}^{t}\frac{Bg^{3/2}(s)-Af^{3/2}(s)}{t^{2}}s^{2}dsdt\right)$
(121)

For $F=\frac{12\pi}{5k_{p}}\left(q^{2}-Gm_{p}^{2}\right)$,
$E=\frac{12\pi}{5k_{p}}\left(q^{2}+Gm_{p}m_{e}\right)$,
$B=\frac{12\pi}{5k_{e}}\left(q^{2}-Gm^{2}_{e}\right)$, and
$A=\frac{12\pi}{5k_{e}}\left(q^{2}+Gm_{p}m_{e}\right)$. It is easy to see that
any such $\hat{f}$ and $\hat{g}$ are $C^{2}$, and that $T$ maps $X$ into
itself if $\delta$ is small enough. For we can say

$\displaystyle|\alpha-\hat{f}(r)|\leq[F(\alpha+\epsilon)^{3/2}+E(\beta+\epsilon)^{3/2}]\int_{0}^{r}\int_{0}^{t}\frac{s^{2}}{t^{2}}dsdt=\frac{F(\alpha+\epsilon)^{3/2}+E(\beta+\epsilon)^{3/2}}{6}r^{2}$
(122)

which can be made smaller than $\epsilon$ for $\delta$ small enough. This is a
crude estimate, but for fixed $\alpha,\beta$, and $\epsilon$, it gives a lower
bound on the length of the interval. Note that, for fixed $\epsilon$, this
lower bound varies continuously in $\beta$ and $\alpha$, so that for small
changes in $(\alpha,\beta)$, we may assume some minimal interval length. All
of this is of course also true for $\hat{g}$. So $T$ maps $X$ to itself. It is
also straightforward to check that a fixed point of (121) is a solution to
(118) and (119). Then we next must show that $T$ is a contraction on $X$. This
can be seen from

$\displaystyle|\hat{f}_{1}(r)-\hat{f}_{2}(r)|$
$\displaystyle=\left|\int_{0}^{r}\int_{0}^{t}\frac{(Ff_{1}^{3/2}(s)-Eg_{1}^{3/2}(s))-(Ff_{2}^{3/2}(s)-Eg_{2}^{3/2}(s))}{t^{2}}s^{2}dsdt\right|$
$\displaystyle\leq\int_{0}^{r}\int_{0}^{t}\frac{|(Ff_{1}^{3/2}(s)-Eg_{1}^{3/2}(s))-(Ff_{2}^{3/2}(s)-Eg_{2}^{3/2}(s))|}{t^{2}}s^{2}dsdt$
$\displaystyle\leq||(f_{1},g_{1})-(f_{2},g_{2})||_{X}(F+E)(\alpha+\beta+2\epsilon)^{1/2}\frac{r^{2}}{4}$

which is a contraction for $\delta$ small enough. We can say the same thing
about $g$. Then we map apply the Banach Fixed Point Theorem to conclude that
we have a unique classical solution to (118) and (119).

To prove that the solutions depend continuously on the initial data, we need
to show that given some $(\alpha,\beta)$ and $\eta>0$, there is a $\delta$ and
$r_{0}>0$ such that
$||(\alpha,\beta)-(\hat{\alpha},\hat{\beta})||_{\infty}<\delta$ implies
$||(f,g)-(\hat{f},\hat{g})||_{\infty}<\eta$ on $[0,r_{0}]$, where $(f,g)$ and
$(\hat{f},\hat{g})$ are the unique solutions corresponding to the initial data
$(\alpha,\beta)$ and $(\hat{\alpha},\hat{\beta})$, respectively. For any
$(\alpha,\beta)$ both greater than zero, by the above statements, we know that
if $\delta$ is small enough, we can assume some interval $[0,r_{1}]$ on which
they both exist, and
$\max(||f-\alpha||_{\infty},||g-\beta||_{\infty},||\hat{f}-\hat{\alpha}||_{\infty},||\hat{g}-\hat{\beta}||_{\infty})<\epsilon$.
Then we may say that $|f(r)-\hat{f}(r)|<2\epsilon+\delta$, and the same for
$g$. This gives us that

$|f(r)-\hat{f}(r)|<\delta+\frac{r^{2}}{4}(\alpha+\beta+\epsilon+\delta)^{1/2}[F+E](2\epsilon+\delta),$
(123)

which we can make smaller than any $\eta$ by choosing $r_{0}$ and $\delta$
small enough. The same is true for $g$, so we find that the function values,
at least in some nonzero length interval, depend continuously on the initial
data.

Of course, this only holds for a small interval around zero. But once we
consider intervals $[R_{0},R_{1}]$ for $R_{0}>0$, we can use more standard
theorems as now the $\frac{2}{r}$ is bounded below and the left hand side of
the ODE is Lipschitz. Thus we may conclude that in the bulk problem, there is
a unique solution to the ODE system, and the densities at a particular
distance from the origin $u_{e}(r)$ and $u_{p}(r)$ depend on the initial
values in a continuous way. ∎

Thus we have local well-posedness, so we now turn to global existence. For the
solution to a given $(\alpha,\beta)$, define
$\psi(r)=-Au^{3/2}_{p}(r)+Bu_{e}^{3/2}(r)$ and
$\phi(r)=-Eu_{e}^{3/2}(r)+Fu^{3/2}_{p}(r)$ so that (118) and (119) can be
rewritten as

$u^{\prime\prime}_{e}(r)+\frac{2}{r}u^{\prime}_{e}(r)=\psi(r)$ (124)

and

$u^{\prime\prime}_{p}(r)+\frac{2}{r}u^{\prime}_{p}(r)=\phi(r)$ (125)

respectively.

######  Lemma 5.2.

Only solutions such that $\phi(0),\psi(0)<0$ can have nonincreasing densities.

###### Proof.

In this proof, and in many of the subsequent proofs, we use a sort of crude
maximum principle. Suppose that $\phi(0)>0$ and $\psi(0)<0$. Since
$\lim_{r\rightarrow 0}\frac{u^{\prime}_{e}(r)}{r}=\frac{2}{3}\psi(0)$, and
$\lim_{r\rightarrow 0}\frac{u^{\prime}_{p}(r)}{r}=\frac{2}{3}\phi(0)$ we can
say that $u^{\prime\prime}_{p}(0)>0$ and $u^{\prime\prime}_{e}(0)<0$. Then by
continuity, $u_{p}$ is increasing and $u_{e}$ is decreasing on a short
interval. But then on the same interval, $\phi$ is increasing and $\psi$ is
decreasing. By our local existence result, let us continue this solution until
either $u_{e}^{\prime}(r)=0$ or $u_{p}^{\prime}(r)=0$; the conclusion will be
the same if neither happens. In the former case, we can see from (124) that
$u^{\prime\prime}_{e}(r)<0$, so $u_{e}^{\prime}$ will continue decreasing,
while (125) tells us that in the latter case $u_{p}$ will continue increasing.
So such starting values give solutions for which $u_{e}$ is nonincreasing and
$u_{p}$ is nondecreasing, at least while both are positive. We have a similar
outcome if $\phi(0)<0$ and $\psi(0)>0$.

To get $\phi(0)>0$ and $\psi(0)>0$, we must have

$\psi(0)=-A\alpha^{3/2}+B\beta^{3/2}>0$ (126)

and

$\phi(0)=-E\beta^{3/2}+F\alpha^{3/2}>0$ (127)

Combining these inequalities leads to the inequality $BF>EA$, but these are
all physical constants, and it can be checked that this inequality is not
satisfied for our problem. Therefore, we may discount this case.

There are some boundary cases we also need to consider. The case
$\phi(0)=\psi(0)=0$ and the cases $\psi(0)=0$, $\phi(0)>0$ and $\psi(0)>0$,
$\phi(0)=0$ all cannot happen for the same reason that $\psi(0)>0$,
$\phi(0)>0$ cannot happen. The case $\phi(0)<0$, $\psi(0)=0$ becomes the same
as the case $\phi(0)<0$, $\psi(0)>0$, while the case $\phi(0)=0$, $\psi(0)<0$
becomes the same as the case $\phi(0)>0$, $\psi(0)<0$. ∎

Although we can find solutions to (118) and (119) with initial conditions such
that either $\phi(0)$ or $\psi(0)$ is not negative (the situation considered
in the previous lemma), we will show below in the subsection on the
atmospheric problem that such solutions cannot be integrable.

We are left to consider cases which satisfy

$\phi(0)=-E\beta^{3/2}+F\alpha^{3/2}<0$ (128)

and

$\psi(0)=-A\alpha^{3/2}+B\beta^{3/2}<0$ (129)

For this to be satisfied, we need to have $EA>BF$, which we know is true.
Similar to the ratio of the total number electrons to protons, (128) and (129)
severely restrict the cases under consideration, as we must have both
$\alpha<(E/F)^{2/3}\beta$ and $\alpha>(B/A)^{2/3}\beta$. In fact, if we write
these out, we see they are essentially (LABEL:bounds):

$(q^{2})^{2/3}\left[\frac{1-\frac{Gm_{q}^{2}}{q^{2}}}{1+\frac{Gm_{e}m_{p}}{q^{2}}}\right]^{2/3}<\frac{\alpha}{\beta}<(q^{2})^{2/3}\left[\frac{1+\frac{Gm_{e}m_{p}}{q^{2}}}{1-\frac{Gm_{p}^{2}}{q^{2}}}\right]^{2/3}$
(130)

In describing the global structure of the solutions to the bulk problem, we
may ask whether there is a “special solution”, similar to the “no atmosphere”
solution presented in [9]: that is, a solution to this system such that the
densities are proportional to each other. In fact, we can construct such
solutions in the general case where we have the system

$\Delta(u^{(d-3)/3})=-Av+Bu$ (131)

and

$\Delta(v^{(d-3)/3})=-Eu+Fv$ (132)

for $d>3$ (still using the physical constants).

######  Lemma 5.3.

For $d>3$, the system given by (131), (132), and

$u(0)=\alpha>0,v(0)=\beta>0,\text{ and }u^{\prime}(0)=v^{\prime}(0)=0;$ (133)

has a solution such that $u$ is proportional to $v$. If $d>\frac{18}{5}$,
there is such a solution with compact support. If $3<d<6$, there is only one
such solution.

###### Proof.

Note that this is the system of PDEs we would get from the Euler-Lagrange
equations of our energy functional if our kinetic energy term was proportional
to $\int_{\mathbb{R}^{3}}u^{d/3}$. Let us assume $u=k_{d}v$. Plugging this
into (131) and (132), we get a system which allows us to solve for $k$. We
find that our desired $k_{d}$ is a solution to

$H_{d}(k)=k^{d/3}F+Ak-Ek^{(d-3)/3}-B=0$ (134)

For $3<d<6$, $H_{d}$ has a single root, which is easy to see as $H_{d}(0)<0$
while $\lim_{k\rightarrow\infty}H_{d}(k)=\infty$ and $H_{d}^{\prime\prime}>0$.

Now that we have a $k_{d}$ for each $d$, we need to solve $\Delta
v^{(d-3)/3}=(-E+k_{d}F)v$, where we will assume spherical symmetry. As above,
let us set $v^{(d-3)/3}=\nu$. Then in spherical coordinates, this equation
becomes

$\nu^{\prime\prime}(r)+\frac{2}{r}\nu(r)=(-E+k_{d}F)\nu^{3/(d-3)}$ (135)

How we can treat this equation depends on the sign of $-E+k_{d}F$. We can
determine this by inserting $E/F$ into $H_{d}(k)$. We get

$\left(\frac{E}{F}\right)^{d/3}F+A\frac{E}{F}-E\left(\frac{E}{F}\right)^{(d-3)/3}-B=A\frac{E}{F}-B>0.$

Since $H_{d}$ is convex and has a single root, we can then conclude that
$E/F>k_{d}$ which implies that $-E+k_{d}F<0$. But this means that (135) can be
transformed into the Lane-Emden equation of index $3/(d-3)$. For the range
$3<d\leq 6$, this covers indices of $[1,\infty)$. It is well known that the
Lane-Emden equation has solutions, although only for indices less than or
equal to 5 are they integrable. ∎

The Lane-Emden equation has analytic solutions when the index is 1 or 5 (also
0, but this is not in our range), which correspond to $d=6$ or
$d=\frac{18}{5}$. The case under consideration here is for $d=5$, while [9]
presents the linear case $d=6$. As has been mentioned, [9] discusses $d=6$
fully. The $d=\frac{18}{5}$ case has found some applications in General
Relativity, but we will not consider it further. Our $d=5$ case corresponds to
an index of $3/2$, solutions of which are known to have compact support.

As the physical interpretation of the Lane-Emden equation is the density of a
self-gravitating polytropic fluid, we can then interpret our solutions which
are close to the special solutions as perturbations of this. As a follow up to
the zero-G model discussed in section 2, one can easily see what happens to
the polytropic solutions when $G\rightarrow 0$. Setting $G=0$ means that $E=F$
and $B=A$. As $H_{d}(k)$ has a single positive root, we must have $k_{d}=1$.
So as $G\rightarrow 0$, $k_{d}\rightarrow 1$. Therefore, when we look for the
special solutions of the zero-G model, we get $u=v$, and the electric energy
also vanishes. This leaves only the kinetic energy. Obviously, this has no
nonzero minimizer, so we get the same result as above.

Since we would expect the proton density to be almost equal to the electron
density, we may seek solutions of the form $u=v$ as approximations. These will
not be solutions of (131) and (132); even in the case when $G=0$, this
produces $u=v=\gamma$ for $\gamma$ any constant, which clearly is not
integrable. If we do want densities of this form, we instead minimize the
energy functional

$E_{n}(u)=(k_{e}+k_{p})\int_{\mathbb{R}^{3}}u^{d/3}(s)d^{3}s-G(m_{e}+m_{p})^{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{u(s)u(t)}{|s-t|}d^{3}sd^{3}t$
(136)

The Euler-Lagrange equation is

$-\frac{d}{3}(k_{e}+k_{p})u^{(d-3)/3}(s)-G(m_{e}+m_{p})^{2}\int_{\mathbb{R}^{3}}\frac{u(t)}{|s-t|}d^{3}t+\lambda=0$
(137)

So applying $-4\pi\Delta$ as above, and letting $u^{(d-3)/3}=\nu$, we get

$\frac{d}{3}(k_{e}+k_{p})\Delta\nu(s)=-4\pi
G(m_{e}+m_{p})^{2}\nu^{3/(d-3)}(s)$ (138)

But if we assume spherical symmetry, this again gives us a (scaled) Lane-Emden
equation of index $(d-3)/3$. This scaling will be different from that for the
special solution of (131) and (132), so the densities will be different
although have the same shape. This locally neutral approximation is precisely
what has been done for decades, and is exactly the assumption that we are
removing in our model.

Returning to our main problem, from our assumption that $\phi(0)<0$ and
$\psi(0)<0$, we see that both $u_{e}$ and $u_{p}$ in a local neighborhood of 0
are nonincreasing. Applying our local existence result, we continue solving
outwards until at least one of $u_{p}$ or $u_{e}$ vanishes or has positive
slope. Before we consider what cases this gives us, we must show that these
are the only possibilities.

######  Lemma 5.4.

There is no solution such that $\phi(0)<0$, $\psi(0)<0$, and both densities
are nonincreasing with unbounded support.

###### Proof.

Let $\lim u_{p}=u_{p}^{\infty}$ and $\lim u_{e}=u_{e}^{\infty}$. From (124)
and (125) we see that both $\psi$ and $\phi$ must have a limit of zero. This
is only possible if both $u_{p}^{\infty}$ and $u_{e}^{\infty}$ are zero. To
eliminate this possibility, we take advantage of the special solutions
discussed above. Suppose that $(v_{1},u_{1})$ is a solution of this form. Then
we can find a special solution $(v_{2},u_{2})$ such that either
$u_{1}(0)=u_{2}(0)$ and $v_{1}(0)\neq v_{2}(0)$ or $u_{1}(0)\neq u_{2}(0)$ and
$v_{1}(0)=v_{2}(0)$; uniqueness prevents both equalities.

There are then four cases, but all can be treated in the same way. Consider
the case that $u_{1}(0)=u_{2}(0)$ and $v_{1}(0)>v_{2}(0)$. We compute

$(u_{1}-u_{2})^{\prime\prime}(0)=F(u_{1}^{3/2}(0)-u_{2}^{3/2}(0))-E(v_{1}^{3/2}(0)-v_{2}^{3/2}(0))<0$
(139)

so that $(u_{1}-u_{2})$ is decreasing for a short interval around the origin,
implying that on that interval, $u_{2}>u_{1}$. We also compute

$(v_{1}-v_{2})^{\prime\prime}(0)=B(v_{1}^{3/2}(0)-v_{2}^{3/2}(0))-A(u_{1}^{3/2}(0)-u_{2}^{3/2}(0))>0$
(140)

So $v_{1}$ is initially greater than $v_{2}$, and this difference is on some
interval increasing. Now we apply the existence result to continue these
solutions until either $u_{1}^{\prime}=u_{2}^{\prime}$ or
$v_{1}^{\prime}=v_{2}^{\prime}$; it is clear that $u_{1}<u_{2}$ and
$v_{1}>v_{2}$ at any distance less than or equal to this. If this never
happens, we get a contradiction since $u_{1}<u_{2}$, but $u_{2}$ vanishes at
some finite distance. If $u_{1}^{\prime}=u_{2}^{\prime}$, (119) gives us that
$(u_{1}-u_{2})^{\prime\prime}<0$, and if $v_{1}^{\prime}=v_{2}^{\prime}$,
(118) gives us that $(v_{1}-v_{2})^{\prime\prime}>0$. So $u_{1}\leq u_{2}$,
which again gives a contradiction since $u_{2}$ has compact support. The other
three cases are treated similarly, and allow us to conclude that there is no
solution to the bulk problem with unbounded support. ∎

This same technique can be used to show that the special solutions are the
only solutions such that the densities vanish at the same point.

######  Lemma 5.5.

For $d=5$, the special solution given by lemma 5.3 is the only solution such
that the densities have the same support.

###### Proof.

Suppose $(u_{1},v_{1})$ is another such solution. We can again find
$(u_{2},v_{2})$ such that either $u_{1}(0)=u_{2}(0)$ and $v_{1}(0)\neq
v_{2}(0)$ or $u_{1}(0)\neq u_{2}(0)$ and $v_{1}(0)=v_{2}(0)$. Assuming the
case that $u_{1}(0)=u_{2}(0)$ and $v_{1}(0)>v_{2}(0)$, we can again conclude
that $u_{1}\leq u_{2}$ and $v_{1}\geq v_{2}$. In fact, since we can say that
$u_{1}^{\prime}\leq u_{2}^{\prime}$ and at some point
$u_{1}^{\prime}<u_{2}^{\prime}$, we can conclude that $u_{1}<u_{2}$. With the
same reasoning $v_{1}>v_{2}$. But since $v_{2}$ and $u_{2}$ vanish at the same
point, we have a contradiction. As the other cases are similar, we conclude
that the special solutions are the unique solutions such that the densities
vanish at the same point. ∎

Then we have returned to our stopping criteria. In the simplest case, one of
the densities vanishes first, and we find ourselves in the exterior problem,
to be elaborated on below. The case when both vanish at the same distance must
be a special solution, as shown above.

######  Lemma 5.6.

If $u_{f}^{\prime}(R_{0})>0$ at $R_{0}>0$, then for $R>R_{0}$,
$u_{f}^{\prime}\geq 0$ while $u_{\sim f}^{\prime}\leq 0$.

###### Proof.

Let us assume that $u_{p}^{\prime}$ becomes positive before $u_{e}^{\prime}$.
By continuity, $u_{p}^{\prime}$ must first be zero, so let us consider what
can happen at $R_{0}$, the last point at which $u_{p}^{\prime}=0$ before it
becomes positive. We are going to examine all the possibilities for different
values of $u_{p}$, $u_{e}$, $u_{p}^{\prime}$, and $u_{e}^{\prime}$ at $R_{0}$
and show that all have the same behavior.

First, consider if $u_{e}^{\prime}(R_{0})=u_{p}^{\prime}(R_{0})=0$. We must
always have one of $\psi$ or $\phi$ be strictly negative, and since we are
assuming $u_{p}^{\prime}$ becomes positive first, this must be $\psi$.
Therefore $u_{e}^{\prime\prime}(R_{0})<0$, and on some interval
$[R_{0},R_{0}+\delta]$, $u_{e}^{\prime}$ is negative. If $\phi(R_{0})>0$, then
$u_{p}^{\prime}$ is positive on an interval also. If $\phi(R_{0})=0$ then
$u_{p}^{\prime\prime}(R_{0})=0$, and we use $\psi(R_{0})<0$ to conclude that
$u^{\prime\prime}(R_{0})<0$. But this means that
$\phi^{\prime\prime}(R_{0})>0$. Since $\phi^{\prime}(R_{0})=0$, we again
conclude that $u_{p}$ is increasing while $u_{e}$ is decreasing on some
interval $[R_{0},R_{0}+\delta]$.

Now we may consider when $u_{p}^{\prime}(R_{0})=0$ and
$u^{\prime}_{e}(R_{0})<0$. If $\phi(R_{0})>0$ then
$u_{p}^{\prime\prime}(R_{0})>0$ and again we find that on some interval
$[R_{0},R_{0}+\delta]$ $u^{\prime}_{e}$ is decreasing while $u^{\prime}_{p}$
is increasing. If $\phi(R_{0})=0$, then $u_{p}^{\prime\prime}=0$. We must then
account for the behavior of $u_{e}$. Since $u_{e}^{\prime}(R_{0})<0$, then
$\phi^{\prime}(R_{0})>0$. This implies that $u_{p}^{\prime}\geq 0$ on an
interval, and we come out with the same result.

So we can conclude that their is some interval $[R_{0},R_{0}+\delta]$ such
that $u_{p}$ is increasing and $u_{e}$ is decreasing. It remains to conclude
that, as long as (118) and (119) are valid (when both densities are positive),
$u_{p}$ is nondecreasing and $u_{e}$ is nonincreasing. Note that in all cases
above, on this same short interval, $\phi\geq 0$ and $\psi\leq 0$.
$u_{p}^{\prime}$ cannot become negative without first being zero and
$u_{e}^{\prime}$ cannot become positive without first being zero. Suppose we
continue our solutions until the first time one of these happens, again
suppose it is $u_{p}^{\prime}(R_{1})=0$. Then $u_{p}^{\prime\prime}(R_{1})$
will have the same sign as $\phi$. But $\phi$ cannot decrease since $u_{e}$ is
decreasing and $u_{p}$ is increasing on the interval $[R_{0},R_{1}]$.
Therefore $\phi(R_{1})\geq 0$, and we can conclude that $u_{p}^{\prime}$
cannot become negative. The same argument is true for $u_{e}^{\prime}$ not
becoming positive. Therefore, for $R>R_{0}$, $u_{e}^{\prime}(R)\leq 0$ and
$u_{p}^{\prime}(R)\geq 0$. Obviously, the same result holds if
$u_{e}^{\prime}$ first becomes positive. ∎

To conclude our classification of the bulk problem, we must then determine
what happens when either $u_{p}^{\prime}$ or $u_{e}^{\prime}$ becomes
positive.

######  Lemma 5.7.

If either density is at some point increasing, the other density will vanish,
and the increasing density will remain bounded on the vanishing density’s
support.

###### Proof.

We have already seen one density becomes nonincreasing while the other is
nondecreasing. Only three things can happen from here. For simplicity, assume
that $u_{p}^{\prime}>0$. In the first situation, $u_{p}$ stays bounded while
$u_{e}$ vanishes with compact support; in the second situation, $u_{p}$
develops a singularity while $u_{e}>0$; and in the third situation, both
$u_{p}$ and $u_{e}$ have support on the entire line. We will here prove that
the second and third options cannot happen. To show that the second situation
cannot happen, we multiply (118) by $A/F$ and add this to (119). This gives us

$(u_{e}^{\prime\prime}+\frac{A}{F}u_{p}^{\prime\prime})+\frac{2}{r}(u_{e}^{\prime}+\frac{A}{F}u_{p}^{\prime})=(B-\frac{A}{F}E)u_{e}^{3/2}$
(141)

If $u_{p}$ developed a singularity at some $R_{0}$, then the left hand side
would become unbounded as $r\rightarrow R_{0}$, while the right hand side
would remain bounded, a contradiction. To see that the third option cannot
happen, first note that this case implies that $\psi$ has a limit of zero. We
can therefore conclude that $u_{p}$ must be bounded and has some finite limit
as $r\rightarrow\infty$. This implies that
$\lim_{r\rightarrow\infty}\phi(r)=(\frac{B}{A}F-E)\lim_{r\rightarrow\infty}u_{e}^{3/2}(r)$.
Since $(\frac{B}{A}F-E)<0$ but $\phi\geq 0$, we must have
$\lim_{r\rightarrow\infty}u_{e}^{3/2}(r)=0$. Of course, this implies
$\lim_{r\rightarrow\infty}u_{p}^{3/2}(r)=0$, a contradiction. Then we have
eliminated the second and third scenarios, so only the first scenario is
possible. ∎

Let us summarize our results for the bulk problem:

  1. 1.

We cannot have both densities increasing. This follows from lemma 5.2 and
lemma 5.6

  2. 2.

If a density is increasing, it is bounded on the decreasing density’s support,
which must be compact. This is lemma 5.7

  3. 3.

The only solutions for which the densities have the same support are the
special solutions. This is lemma 5.5

So given any initial values for the bulk problem, if we do not have a special
solution, we see that one of the densities will vanish at some finite radius
$R_{0}$. This transfers directly to the atmospheric problem.

###  5.2 The Atmospheric Problem

For the atmospheric problem, we are trying to solve the initial value problem

$\displaystyle\Delta u$ $\displaystyle=Du^{3/2}\text{ on }(R_{0},R_{1})$ (142)
$\displaystyle u(R_{0})$ $\displaystyle=\gamma>0$ (143)

Where $D>0$, $R_{0}>0$, and $R_{1}$ is either where $u$ vanishes, or $\infty$.
With radial symmetry, it is a standard exercise to show that this problem is
locally well-posed. We therefore focus only on the global structure of the
solutions.

There are a few simple observations which can be made. First, if $u^{\prime}$
becomes zero on the support of $u$, then $u$ will always be nondecreasing
after that point. This is easy to see as if $u^{\prime\prime}=Du^{3/2}>0$, the
first derivative must be increasing. This tells us that we cannot have any
oscillatory behavior and that any solution which ever hits zero must have
always been decreasing. When we combine the bulk and atmospheric problem, this
will also prove the statement that we can ignore any solution of the bulk
problem for which one of the densities is nondecreasing.

Second, let us consider what the limiting values of $u$ may be. In particular,
we want to conclude that there are no solutions such that
$u(r)\rightarrow\alpha>0$ as $r\rightarrow\infty$ (since there is no
oscillatory behavior, we know $u$ is monotonic, so such a limit exists). For
suppose that such a solution exists, then $u$ may approach $\alpha$ from above
or below. In either case, we must have $u^{\prime}(r)\rightarrow 0$, so for
$r$ larger than some $R$, we can conclude that
$u^{\prime\prime}(r)>\alpha-\epsilon$>0 contradicting that
$u^{\prime}(r)\rightarrow 0$. So if $u$ remains positive, it can only have
limits of $\infty$ or 0.

Further, there may be only one strictly positive solution which approaches
zero at infinity. For suppose that $u_{1}$ and $u_{2}$ were two such
solutions, and let $u_{1}^{\prime}(R_{0})>u_{2}^{\prime}(R_{0})$. Then
$u_{1}>u_{2}$ on $(R_{0},R]$ for any finite $R$ by the maximum principle, and
the difference is largest at $R$. But by assumption, if we take $R$ large
enough we must have both $u_{1}(r)$ and $u_{2}(r)$ in $(0,\epsilon)$ for any
$\epsilon>0$. Thus we can force $u_{1}^{\prime}(R_{0})$ to be arbitrarily
close to $u_{2}^{\prime}(R_{0})$, so we must have equality and $u_{1}=u_{2}$.
Call this possible solution $u_{\infty}$. This shows that for fixed $R_{0}$
and $u(R_{0})$, if $u$ is a solution such that
$u^{\prime}(R_{0})>u^{\prime}_{\infty}(R_{0})$, $u$ must be unbounded.

Finally, concerning solutions which have compact support, we can say that if
we have two solutions $u_{1}$ and $u_{2}$ satisfying
$u_{1}(R_{0})=u_{2}(R_{0})$, the one with the larger slope at $R_{0}$ must hit
zero at a larger $r$. This follows from a maximum principle as the solution
with a larger initial slope must have a larger function value than the other
solution, at least as long as one of the functions is positive. This tells us
that the initial slopes have a linear ordering: larger initial slopes lead to
larger radii when the density becomes zero, if it becomes zero at all.
However, nothing above tells us that any solution ever hits zero, so that is
what we turn to now.

If we pose this question as the boundary value problem,

$\displaystyle\Delta u$ $\displaystyle=f(u)\text{ on }\Omega$ (144)
$\displaystyle u$ $\displaystyle=g\text{ on }\partial\Omega$ (145)

where $\Omega$ is the annular domain with inner radius $R_{0}$ and outer
radius $R$, $g=a>0$ on the inner radius and $g=0$ on the outer radius, it is a
standard exercise that a classical unique solution exists. Since the boundary
values are spherically symmetric, by uniqueness, so is the solution.
Therefore, the solution found here will be the same as the unique solution to
the ODE problem with the corresponding $u^{\prime}(R_{0})$ value.

From the above arguments, we may conclude that given a fixed $a>0$ and
$R_{0}>0$, for every $R>R_{0}$, we can find a unique strictly decreasing
solution $u_{R}$, and that $u_{R_{1}}(r)\geq u_{R_{2}}(r)$ if $R_{1}>R_{2}$,
with the inequality strict for $R_{0}<r\leq R_{2}$. Now, we want to prove that
there is a solution (we have already stated it must be unique) with unbounded
support which has $\lim_{r\rightarrow\infty}u(r)\rightarrow 0$.

######  Theorem 5.8.

For any given $R_{0}$ and $a>0$, there is a unique spherically symmetric
solution to $\Delta u=Du^{3/2}$ which is both strictly positive and satisfies
$\lim_{r\rightarrow\infty}u(r)=0$.

###### Proof.

First, since each $u_{R}$ is strictly decreasing and $u_{R_{1}}(r)\geq
u_{R_{2}}(r)$ if $R_{1}>R_{2}$, $\lim_{R\rightarrow\infty}u_{R}(r)$ exists for
every $r$. Define this value to be $\hat{u}(r)$. The claim is that
$\hat{u}(r)$ is the desired function. To show this, we are going to show that
$u_{R}\rightarrow\hat{u}$ in $C^{2}{[R_{0},R_{1}]}$ for any
$R_{0}<R_{1}<\infty$. Then in particular we will have that
$D^{\alpha}u_{R}(r)\rightarrow D^{\alpha}\hat{u}(r)$ for any $|\alpha|\leq 2$
and any $R_{0}<r<\infty$. So $\hat{u}$ will solve the ODE.

To get this convergence, we use Rellich-Kondrachov to get a compact embedding
of the set $\\{u_{R}\\}_{R>N}$ when we have restricted to
$\Omega_{N}=(R_{0},N)$. So fix $N$ and consider the set $\\{u_{R}\\}_{R>N}$.
In particular, each of these functions is strictly greater than zero on
$\Omega_{N}$. Then applying some elliptic regularity estimates, we conclude
that each $u_{R}$ is in $H^{3}(\Omega_{N})$ (meaning the $L^{2}$ Sobolev
space) and further satisfies

$||u_{R}||_{3,\Omega_{N}}\leq
C(||u_{R}^{3/2}||_{1,\Omega_{N}}+||u_{R}||_{0,\Omega_{N}}+||u_{R}||_{3,\partial\Omega_{N}})$
(146)

Thinking of $\Omega_{N}$ as a set in $\mathbb{R}^{3}$, the last term is
bounded as $u_{R}$ is less or equal to $a$ on the boundary. The middle term is
also bounded by $a$ multiplied by the size of $\Omega_{N}$. To bound the first
term, we need to have some control over $u_{R}^{\prime}$, and to do this, we
apply the same estimate one derivative lower:

$||u_{R}||_{1,\Omega_{N}}\leq||u_{R}||_{2,\Omega_{N}}\leq
C(||u_{R}^{3/2}||_{0,\Omega_{N}}+||u_{R}||_{0,\Omega_{N}}+||u||_{2,\partial\Omega_{N}})\leq
K$

for some fixed $K$, for all $R>N$. This allows us to get a uniform bound on
(146). So we not only have that $\\{u_{R}\\}_{R>N}$ is contained in
$H^{3}(\Omega_{N})$, we have that it is bounded as well. Then, treating these
as functions on the real line, Rellich-Kondrachov gives us that this set
embeds compactly into $C^{2}(\overline{\Omega})=C^{2}([R_{0},N])$. In
particular we get pointwise convergence as described above, so $\hat{u}$ is a
solution to the ODE. ∎

$\hat{u}$ is then a sort of boundary case to solutions of (142): solutions
with $u^{\prime}(R_{0})<\hat{u}^{\prime}(R_{0})$ have compact support while
solutions with $u^{\prime}(R_{0})>\hat{u}^{\prime}(R_{0})$ are unbounded.
Although not interesting from a physical perspective, we can say more about
these unbounded solutions.

######  Lemma 5.9.

Any solution greater than $\hat{u}$ develops a singularity.

###### Proof.

The key here is what we have already proven about these solutions: they cannot
remain bounded as $r\rightarrow\infty$. Then given any such solution of this
type, we may assume that for large enough $r$, $u(r)>\gamma$ for some large
$\gamma$. From the ODE, we may also conclude that $u^{\prime}$ cannot be
bounded. Fix a large $R$, and consider $f_{R}(r)=\frac{\omega}{(R+1-r)^{4}}$.
If $\omega$ is chosen large enough, then on the interval $[R,R+1)$,

$\Delta f_{R}(r)=\frac{20\omega}{(R+1-r)^{6}}+\frac{8\omega}{r(R+1-r)^{5}}\leq
D\frac{\omega^{3/2}}{(R+1-r)^{6}}=Df_{R}^{3/2}(r)$ (147)

Note that $\omega$ is chosen independently of $R$. We have
$\Delta(u(r)-f_{R}(r))=D(u^{3/2}{r}-f_{R}^{3/2}(r))>0$ on a small interval
$[R,R+\delta]$ if $u(R)>f_{R}(R)$, and $u^{\prime}(R)>f_{R}^{\prime}(R)$. This
implies that $u(r)>f_{R}(r)$ as long as (147) holds. But since $u(r)$ and
$u^{\prime}(r)$ are unbounded, we may find an $R$ so that $u(R)>f_{R}(R)$ and
$u^{\prime}(R)>f^{\prime}_{R}(R)$. Therefore $u(r)>f_{R}(r)$ on $[R,R+1)$, so
$u(r)$ must develop a singularity. ∎

A final point we want to make is to show that $\hat{u}$ is integrable.

######  Lemma 5.10.

$\hat{u}$ is integrable.

###### Proof.

Recall that for each $R>0$ and each $a>0$ we get such a $\hat{u}$, so now we
label them as $\hat{u}_{R,a}$. We will use that we have a special solution:
$v(r)=\frac{C}{r^{4}}$ for $C=144/D^{2}$. Obviously, since the equation is not
linear, we cannot scale this to get any other solutions. But we see that if
$v(R)=a$, then $\hat{u}_{R,a}$ coincides with $v(r)$; this must be the case
since we have already shown above that the solution to the ODE which has
unbounded support and converges to zero is unique. For simplicity, fix $R$.
Then we have two scenarios: either $a>v(R)$ or $a<v(R)$.

In the latter case, suppose at some point $R_{1}$
$\hat{u}_{R,a}(R_{1})=v(R_{1})$. By uniqueness of the solution converging to
zero, both of these functions must after $R_{1}$ coincide. Therefore
$\hat{u}_{R,a}\leq v$, and decays at least as fast as $r^{-4}$, which is
integrable outside of a ball on $\mathbb{R}^{3}$. In the former case, we
consider the function $\overline{v}(r)=\frac{\overline{c}}{r^{4}}$, where
$\overline{c}$ is chosen so that $\overline{v}(R)>\hat{u}_{R,a}(R)=a$. Since
$\overline{c}>144/D^{2}$, one can check that, treating the function as
spherically symmetric on $\mathbb{R}^{3}$, $\Delta v\leq D\overline{v}^{3/2}$.
Now suppose that there is some $R_{1}$ so that
$\hat{u}_{R,a}(R_{1})=\overline{v}(R_{1})$. Then, for at least some small
interval $[R_{1},R_{1}+\delta]$, $\hat{u}_{R,a}>\overline{v}$ (we could not
have equality since $\overline{v}$ does not satisfy the ODE). But then we have
on that interval $\Delta(\hat{u}_{R,a}-\overline{v})\geq
D(\hat{u}_{R,a}^{3/2}-\overline{v}^{3/2})>0$. By the maximum principle, this
implies that $\hat{u}_{R,a}(R_{1}+\delta)>\overline{v}(R_{1}+\delta)$, so we
may extend the interval longer and the functions would always get further
apart. But since they both converge to zero, this cannot happen. So we see
that $\hat{u}_{R,a}<\overline{v}$, and therefore goes to zero at least as fast
as $r^{-4}$ and is therefore integrable. So we have shown that in all
scenarios $\hat{u}_{R,a}$ is an integrable function. ∎

###  5.3 Combining the Problems

Now that we have addressed the atmospheric problem and the bulk problem in
isolation, we can bring them together to solve the system of equations. The
method is clear. Given any pair of central densities which satisfy (130), we
use the results of subsection 5.1 to find the unique solution of the system
continuing outward radially until one of the densities vanishes, say at
$R_{0}$. The nonvanishing density will satisfy $u(R_{0})=a>0$ and
$u^{\prime}(R_{0})=b$. We then use these values as our initial values in the
atmospheric problem and apply the results of subsection 5.2.

One issue which needs to be addressed is regularity. We have established that
the solutions of both the bulk and atmospheric problems are $C^{2}$ on their
domains. But the density which is nonzero in the atmosphere might have a jump
in its second derivative at the interface of the two problems. To show that
this is not the case, let us only consider those initial values such that
$u^{3/2}_{p}=\rho_{p}$ and $u^{3/2}_{e}=\rho_{e}$ are integrable, recalling
that in section 2 we transformed the actual densities by the inverse of this.
It is clear that these will be the densities which have compact support, as
well as the limiting case discussed in theorem 5.8. Of course, for our
minimization problem, all the other solutions are nonviable since one of the
densities is unbounded.

Before we proceed, we need the following extension of lemma 4.2 in [8].

######  Lemma 5.11.

Let $f$ be bounded, locally Hölder continuous with exponent $\alpha\leq 1$,
and in $L^{p}$ for any $1\leq p<\infty$ on $\mathbb{R}^{n}$. Then if $w$ is
the Newtonian potential of $f$, $w\in C^{2}(\mathbb{R}^{n})$.

###### Proof.

Consider $f_{R}(y)=\chi_{B_{R}}f(y)$, and call $w_{R}$ its Newtonian
potential. Then $f_{R}$ satisfies the hypotheses of lemma 4.2 in [8] in
$B_{R}$, and so $w_{R}$ is $C^{2}(B_{R})$. By that lemma, we have for any
$x\in B_{R}$

$D_{ij}w_{R}(x)=\int_{B_{R}}D_{ij}\Gamma(x-y)(f_{R}(y)-f_{R}(x))dy-
f_{R}(x)\int_{\partial B_{R}}D_{i}\Gamma(x-y)\nu_{j}(y)ds_{y}$ (148)

where $\Gamma$ is the fundamental solution of Laplace’s equation. For $x$ such
that $|x|<R-1$, rewrite this as

$\displaystyle D_{ij}w_{R}(x)$
$\displaystyle=\int_{B_{1}(x)}D_{ij}\Gamma(x-y)(f_{R}(y)-f_{R}(x))dy+\int_{B_{R}\backslash
B_{1}(x)}D_{ij}\Gamma(x-y)(f_{R}(y)-f_{R}(x))dy$ (149) $\displaystyle-
f_{R}(x)\int_{\partial B_{R}}D_{i}\Gamma(x-y)\nu_{j}(y)ds_{y}$
$\displaystyle=\int_{B_{1}(x)}D_{ij}\Gamma(x-y)(f_{R}(y)-f_{R}(x))dy+\int_{B_{R}\backslash
B_{1}(x)}D_{ij}\Gamma(x-y)f_{R}(y)$ (150) $\displaystyle-
f_{R}(x)\int_{\partial B_{1}}D_{i}\Gamma(x-y)\nu_{j}(y)ds_{y}$

For any selected $x$, we can define $D_{ij}w_{R}(x)$ like this for all
$R>|x|+1$. Now define for $x\in\mathbb{R}^{3}$

$u(x)=\int_{B_{1}(x)}D_{ij}\Gamma(x-y)(f(y)-f(x))dy+\int_{\mathbb{R}^{n}\backslash
B_{1}(x)}D_{ij}\Gamma(x-y)f(y)dy-f(x)\int_{\partial
B_{1}}D_{i}\Gamma(x-y)\nu_{j}(y)ds_{y}$ (151)

This is well defined because we can apply Hölder’s inequality to the second
term. Then for large $R$ we have

$|u(x)-D_{ij}w_{R}(x)|\leq\int_{B^{c}_{R}(x)}\left|D_{ij}\Gamma(x-y)f(y)\right|dy\leq
C\int_{B^{c}_{R}(x)}\frac{f(y)}{|x-y|^{n}}dy$ (152)

which can be made as small as we like if we take $R$ large enough. Then we
have that for any $x\in\mathbb{R}^{n}$, $D_{ij}w_{R}(x)\rightarrow u(x)$.
Further, on any compact set in $\mathbb{R}^{n}$, this convergence is uniform.
Since we also have that $w_{R}(x)\rightarrow w(x)$ uniformly on any compact
set, we can conclude that $u=D_{ij}w$, and that $w$ is
$C^{2}(\mathbb{R}^{n})$. ∎

In the following lemma, we need to use that $\rho_{e}$ and $\rho_{p}$ are
Hölder continuous. For a density with compact support, this will be true as
long as its first derivative does not become unbounded as it approaches the
boundary of its support. This is however ruled out by (118) and (119) since a
large negative first derivative would lead to a positive second derivative.

######  Lemma 5.12.

When $\rho_{e}$ and $\rho_{p}$ are integrable, they are both $C^{2}$ on their
supports.

###### Proof.

When $\rho_{f}$ have compact support, let us define new densities by

$\frac{5}{3}k_{p}\bar{\rho}_{p}^{2/3}(x):=-q^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}-\rho_{e})(y)}{|x-y|}d^{3}y+Gm_{p}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}+m_{e}\rho_{e})(y)}{|x-y|}d^{3}y+\lambda_{p}$
(153)

and

$\frac{5}{3}k_{e}\bar{\rho}_{e}^{2/3}(x):=Cq^{2}\int_{\mathbb{R}^{3}}\frac{(\rho_{p}-\rho_{e})(y)}{|x-y|}d^{3}y+Gm_{e}\int_{\mathbb{R}^{3}}\frac{(m_{p}\rho_{p}+m_{e}\rho_{e})(y)}{|x-y|}d^{3}y+\lambda_{e}$
(154)

where we take (153) to hold on the support of $\rho_{p}$, $\bar{\rho}_{p}$ is
defined to be zero outside of this support, and $\lambda_{p}$ is chosen to
make $\bar{\rho}_{p}$ continuous. $\bar{\rho}_{e}$ is defined analogously. We
can then apply lemma 4.2 in [8] to conclude that both $\bar{\rho}_{e}$ and
$\bar{\rho}_{p}$ are $C^{2}$ on their supports.

Now apply $-\Delta$ to (153) and (154). Suppose that $\rho_{p}$ is the density
which does not vanish at $R_{0}$, but vanishes at $R_{1}$. Considering the
atmospheric and bulk problem separately, we can conclude that

$\Delta(\bar{\rho}_{p}-\rho_{p})=0\text{ on }[0,R_{0}]\text{ and
}[R_{0},R_{1}]$ (155)

Now $\bar{\rho}_{p}(R_{0})=\rho_{p}(R_{0})$, for if not, then we have
$\bar{\rho}_{p}(r)-\rho_{p}(r)=\bar{\rho}_{p}(R_{0})-\rho_{p}(R_{0})$ for
$r<R_{0}$. But neither function has a jump in value at $R_{0}$, so we apply
the maximum principle again in $B_{R_{1}}\backslash B_{R_{0}}$ to conclude
that $\bar{\rho}_{p}(r)-\rho_{p}(r)=\bar{\rho}_{p}(R_{0})-\rho_{p}(R_{0})$ on
this set also. But we know that both functions vanish at $R_{1}$, so we have
that $\bar{\rho}_{p}=\rho_{p}$ on $[0,R_{1}]$. Then $\rho_{p}$ cannot have a
jump in its second derivative as $\bar{\rho}_{p}$ does not.

Using the extension of lemma 4.2 proved above, we can define $\bar{\rho}_{e}$
and $\bar{\rho}_{p}$ in the same way for the case in which one of $\rho_{e}$
or $\rho_{p}$ have noncompact support taking $\lambda_{f}=0$ for the density
with unbounded support, and still conclude that $\bar{\rho}_{e}$ and
$\bar{\rho}_{p}$ are $C^{2}$ on their supports. In the second part of the
argument, assuming $\rho_{p}$ is the density with unbounded support, we have

$\Delta(\bar{\rho}_{p}-\rho_{p})=0\text{ on }[0,R_{0}]\text{ and
}[R_{0},\infty)$ (156)

We can still apply the maximum principle on $[0,R_{0}]$, and we can apply it
to any $[R_{0},R]$ $R<\infty$. So $\bar{\rho}_{p}-\rho_{p}$ is a constant on
$[R_{0},\infty)$. If this constant is not zero then $\bar{\rho}_{p}$ does not
go to zero as $x\rightarrow\infty$, which clearly contradicts its definition.
∎

So the solutions to the system we found are $C^{2}$ on their support, and are
thus solutions to the Euler-Lagrange equations (4) and (5).

##  6 The Structure of Minimizers

We now have quite a few results concerning the minimizers of (3), and here we
combine them to see a more complete picture. First, we have theorems 3.1 and
4.1 which tell us that as long as $N_{e}/N_{p}$ satisfies (25), (3) has a
spherically symmetric minimizer such that the densities integrate to $N_{e}$
and $N_{p}$, and if the inequalities in (25) are sharp, the minimizer has
compact support. Any such minimizer must satisfy the Euler-Lagrange equations,
and we have seen in section 5 that given two central densities
$(\alpha,\beta)$, these equations have a unique solution for which both
densities are decreasing. These results fit together quite well, since we know
that for a fixed $\alpha$, there is a single $\beta$ such that $u_{p}$ has
unbounded support but is integrable, and a single $\beta$ such that $u_{e}$
has unbounded support but is integrable. Such solutions apparently correspond
to the equalities of (25). If we call these two values $\beta_{\alpha}^{h}$
and $\beta_{\alpha}^{l}$, respectively, then every
$\beta\in[\beta_{\alpha}^{l},\beta_{\alpha}^{h}]$ will produce curves with
compact support, and correspond to strict inequalities in (25).

We want to show that the minimizer is unique, and to do this, we must bridge
the gap between the two approaches. It is clear that if we were given a
minimizer for a pair $(N_{e},N_{p})$, we could use the central densities in
solving the Euler-Lagrange equations, and therefore we recover the minimizer.
So to show that the minimizer is unique, we need to show there is a bijection
from $(N_{e},N_{p})$ to $(\alpha,\beta)$.

######  Lemma 6.1.

There is a bijective map from admissible values of $(N_{e},N_{p})$ to
$(\alpha,\beta)$.

###### Proof.

It is clear that given an $(\alpha,\beta)$, we can compute $N_{e}=\int
u_{e}^{3/2}$ and $N_{p}=\int u_{p}^{3/2}$ by solving the ODE system. As a
given pair of central densities has a unique solution, we just need to show
that this map from $(\alpha,\beta)$ to $(N_{e},N_{p})$ is injective. To this
end, we note that solutions to the bulk and atmospheric problem exhibit a
scaling structure. That is, given a solution $(u_{e}(r),u_{p}(r))$, we can
find another solution $(\theta(s),\eta(s))$ such that $\lambda
u_{e}(r)=\theta(s)$ and $\lambda u_{p}(r)=\eta(s)$ where $r=as$ and
$\lambda^{1/2}=a^{2}$.

Then suppose we have $(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$
such that $\alpha_{1}\neq\alpha_{2}$ and $\beta_{1}\neq\beta_{2}$; if just one
of these is an equality, we can use the type of argument given in lemma 5.5 to
show that at least one of $\int u^{3/2}_{e,1}\neq\int u^{3/2}_{e,2}$ or $\int
u^{3/2}_{p,1}\neq\int u^{3/2}_{p,2}$ must be true. Assume that
$\alpha_{1}=\lambda\alpha_{2}$ for $\lambda>1$, so that we have a solution
$(\theta,\eta)$ with $\theta(0)=\alpha_{1}$ and the relation $(\lambda
u_{e,2}(r),\lambda u_{p,2})=(\theta(s),\eta(s))$. Then we have

$\int
u^{3/2}_{e,2}=4\pi\int_{0}^{\infty}r^{2}u^{3/2}_{e,2}(r)dr=4\pi\int_{0}^{\infty}r^{2}u^{3/2}_{e,2}(r)dr=\frac{4\pi
a^{3}}{\lambda^{3/2}}\int_{0}^{\infty}s^{2}\theta^{3/2}(s)ds=a^{-3}\int\theta^{3/2}$
(157)

and similarly,

$\int u^{3/2}_{p,2}=a^{-3}\int\eta^{3/2}.$ (158)

A comparison principle type argument like what was used in lemma 5.5 tells us
that either $\theta\leq u_{e,1}$ and $\eta\geq u_{p,1}$ or $\eta\leq u_{p,1}$
and $\theta\geq u_{e,1}$, depending on the relationship between
$\lambda\beta_{2}$ and $\beta_{1}$. So we can conclude that we cannot have
both $\int u^{3/2}_{e,1}=\int u^{3/2}_{e,2}$ and $\int u^{3/2}_{p,1}=\int
u^{3/2}_{p,2}$. Therefore, the map from the central densities to the number of
particles is injective. ∎

See [9] for further discussion and numerical solutions for these equations.

###  6.1 The Special Relativistic Case

Let us rewrite (15) and (16) as

$\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d}{dr}y_{p}\right)=\frac{1}{\alpha_{p}k_{p}^{3}}\left(1-G\frac{m_{p}^{2}}{q^{2}}\right)(y_{p}^{2}-1)^{3/2}-\frac{1}{\alpha_{p}k_{e}^{3}}\left(1+G\frac{m_{p}m_{e}}{q^{2}}\right)(y_{q}^{2}-1)^{3/2}$
(159)

and

$\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d}{dr}y_{e}\right)=-\frac{1}{\alpha_{e}k_{p}^{3}}\left(1+G\frac{m_{p}m_{e}}{q^{2}}\right)(y_{p}^{2}-1)^{3/2}+\frac{1}{\alpha_{e}k_{e}^{3}}\left(1-G\frac{m_{q}^{2}}{q^{2}}\right)(y_{q}^{2}-1)^{3/2},$
(160)

respectively.

Unlike the proof that minimizers of (3) exist, the proof of their uniqueness
and the relation of the shape of the minimizers to $(N_{e},N_{p})$ does not so
easily generalize to the special relativistic case. The main hindrance in
directly generalizing the results of section 5 is that there are no special
solutions which provide some basis for a type of comparison principle.

######  Lemma 6.2.

There are no decreasing solutions of the form $\rho_{e}=\hat{k}\rho_{p}$.

###### Proof.

This is equivalent to saying there are no solutions of the form
$(y_{p}^{2}-1)=k(y_{q}^{2}-1)$ for
$k=\hat{k}\left(\frac{m_{p}}{m_{e}}\right)^{3}$. So we can conclude
immediately that $k>0$, since we are only interested in positive solutions.
Since we are further interested in decreasing solutions, we can plug
$(y_{p}^{2}-1)=k(y_{q}^{2}-1)$ into the right hand side of (15) to get

$\Delta
y_{p}=\frac{4\pi}{m_{p}c^{2}}\left(\frac{\pi}{3}\right)\left(\frac{2c}{h}\right)^{3}\frac{1}{q^{2}}\left(m_{p}^{3}k\left(1-\frac{Gm^{2}_{p}}{q^{2}}\right)-m_{e}^{3}\left(1+\frac{Gm_{p}m_{e}}{q^{2}}\right)\right)(y_{q}^{2}-1)^{3/2}$
(161)

Since $y_{p}^{\prime}(0)=0$, we can only have a decreasing solution if

$k\leq\frac{m_{e}^{3}}{m_{p}^{3}}\frac{\left(1+\frac{Gm_{p}m_{e}}{q^{2}}\right)}{\left(1-\frac{Gm^{2}_{p}}{q^{2}}\right)}<1$
(162)

So we know that $0<k<1$.

Now we assume such a solution exists and derive a contradiction. Let
$y_{p}=\sqrt{k(y_{q}^{2}-1)+1}$. Then (159) becomes

$k\left[\frac{2}{r}\frac{y_{e}y_{e}^{\prime}}{\sqrt{(y_{q}^{2}-1)k+1}}+\frac{(1-k)(y_{e}^{\prime})^{2}+y_{e}y_{e}^{\prime\prime}(k(y_{q}^{2}-1)+1)}{((y_{q}^{2}-1)k+1)^{3/2}}\right]=(-E+Fk^{3/2})(y_{q}^{2}-1)^{3/2},$
(163)

for $E$ and $F$ appropriately defined. So just as in the proof of lemma (5.3)
we have a compatibility condition on $k$ given by (163) and (160), the latter
of which we rewrite

$y_{e}^{\prime\prime}+\frac{2}{r}y_{e}^{\prime}=(-Ak^{3/2}+B)(y^{2}_{e}-1)^{3/2}$
(164)

This compatibility condition is not as straightforward though. In lemma (5.3),
our condition led us to an algebraic equation we could easily determine had a
solution. Choosing the correct $k$ essentially eliminated one of the
constraint equations, so then we only needed to solve a simple ODE.

To proceed, we note that we can have neither $B=Ak^{3/2}$ nor $E=Fk^{3/2}$,
for in the former case we would have $y_{p}$ be a constant and in the latter
case we would have $y_{e}$ be a constant. Then combine (163) and (164) to get

$\frac{2}{r}\frac{y_{e}y_{e}^{\prime}}{\sqrt{(y_{q}^{2}-1)k+1}}+\frac{(1-k)(y_{e}^{\prime})^{2}+y_{e}y_{e}^{\prime\prime}(k(y_{q}^{2}-1)+1)}{((y_{q}^{2}-1)k+1)^{3/2}}=\frac{-E+Fk^{3/2}}{k(-Ak^{3/2}+B)}\left[y_{e}^{\prime\prime}+\frac{2}{r}y_{e}^{\prime}\right]$
(165)

This can be rearranged to get

$y_{e}^{\prime\prime}+\frac{2}{r}y_{e}^{\prime}=\frac{(y_{e}^{\prime})^{2}}{k(y_{q}^{2}-1)+1}\left[\frac{k(-Ak^{3/2}+B)(1-k)}{((y_{q}^{2}-1)k+1)^{1/2}(-E+k^{3/2}F)-y_{e}k(-Ak^{3/2}+B)}\right].$
(166)

Then we can recombine (166) with (164) to get

$(y_{e}^{\prime})^{2}=\frac{(y_{q}^{2}-1)^{3/2}}{k(1-k)}\left(k(y_{q}^{2}-1)+1\right)\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right].$
(167)

Then we have that

$y_{e}^{\prime}=\frac{(y_{q}^{2}-1)^{3/4}}{\sqrt{k(1-k)}}\left(k(y_{q}^{2}-1)+1\right)^{1/2}\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right]^{1/2}$
(168)

Differentiating (167) we get

$\displaystyle y_{e}^{\prime\prime}$
$\displaystyle=\frac{3(y_{q}^{2}-1)^{1/2}y_{e}}{2k(1-k)}(k(y_{q}^{2}-1)+1)\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right]$
(169)
$\displaystyle+\frac{(y_{q}^{2}-1)^{3/2}y_{e}}{(1-k)}\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right]$
$\displaystyle+\frac{(y_{q}^{2}-1)^{3/2}}{2(1-k)}(k(y_{q}^{2}-1)+1)\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{-1/2}y_{e}-(-Ak^{3/2}+B)\right].$

Then we can combine these last two equations with (164) to get

$\displaystyle-Ak^{3/2}+B$
$\displaystyle=\frac{3y_{e}}{2k(1-k)(y_{q}^{2}-1)}(k(y_{q}^{2}-1)+1)\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right]$
(170)
$\displaystyle+\frac{y_{e}}{(1-k)}\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right]$
$\displaystyle+\frac{k(y_{q}^{2}-1)+1}{2(1-k)}\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{-1/2}y_{e}-(-Ak^{3/2}+B)\right]$
$\displaystyle+\frac{2}{r}\frac{(y_{q}^{2}-1)^{-3/4}}{\sqrt{k(1-k)}}\left(k(y_{q}^{2}-1)+1\right)^{1/2}\left[(-E+Fk^{3/2})(k(y_{q}^{2}-1)+1)^{1/2}-y_{e}k(-Ak^{3/2}+B)\right]^{1/2}$

the right side of which cannot be constant unless $y_{e}$ is constant. ∎

Another major obstacle that prevents a direct generalization of lemma 6.1 is
that, as is easy to tell by direct computation, there is no scaling structure
to these equations comparable to system (118) and (119). Despite these
difficulties, numerics indicate that the structure of the solutions is exactly
the same as in the Newtonian case, once one accounts for the limitations from
(55), see [9].

##  7 Conclusion

In this paper we have studied a two species Thomas-Fermi type model for a non-
rotating non-neutral self-gravitating brown dwarf star in its ground state.
This is an extension of previous results which assumed local neutrality and
therefore were reduced to considering only a one fluid model.

In the Newtonian kinetic energy case, we gave a complete classification of the
unique minimizers of (3) given that the total number of positively and
negatively charged particles satisfied (25). We were also able to show the
existence of minimizers in the special relativistic kinetic energy case.

Of course, these results still make a large number of unrealistic assumptions:
at the very least, actual white dwarfs rotate and are not composed of only
protons and electrons. In future work, we intend to extend this work to models
removing these assumptions. We believe it is however worthwhile to study the
simplest nonneutral case before beginning to address more realistic,
complicated models.

## References

  * [1] Llewellyn Thomas, The calculation of atomic fields, Mathematical Proceedings of the Cambridge Philosophical Society 23(5), (1927).
  * [2] Enrico Fermi, Un Metodo statistico per la determinazione di alcune priorieta dell’atome, Rend. Accad. Naz. Lincei 6, 602–607(1927).
  * [3] Charles Fefferman and Luis Seco, On the energy of a large atom, Bulletin of the American Mathematical Society 23(2), 525–530 (1990).
  * [4] Edward Teller, On the stability of molecules in Thomas-Fermi theory, Rev.Mod. Phys. 34, 627–631 (1962).
  * [5] Lieb, Elliott H., and Walter E. Thirring, Bound for the kinetic energy of fermions which proves the stability of matter, Phy. Rev. Lett. 35, 687–689 (1975).
  * [6] Spruch, Larry, Pedagogic notes on Thomas-Fermi theory (and on some improvements): atoms, stars, and the stability of bulk matter, Rev. Mod. Phys. 63(1) 151–209 (1991).
  * [7] Lieb, Elliott H., and Robert Seiringer, The stability of matter in quantum mechanics. Cambridge University Press, 2010.