Locality of the Thomas-Fermi-von Weizs\"acker Equations

We establish a pointwise stability estimate for the Thomas-Fermi-von Weizs\"acker (TFW) model, which demonstrates that a local perturbation of a nuclear arrangement results also in a local response in the electron density and electrostatic potential. The proof adapts the arguments for existence and uniqueness of solutions to the TFW equations in the thermodynamic limit by Catto et al. (1998). To demonstrate the utility of this combined locality and stability result we derive several consequences, including an exponential convergence rate for the thermodynamic limit, partition of total energy into exponentially localised site energies (and consequently, exponential locality of forces), and generalised and strengthened results on the charge neutrality of local defects.


Introduction
Locality properties of electronic structure models are a key premise in certain state of the art numerical algorithms. A well-established example is "near-sightedness", a locality property of the density matrix which gives rise to linear scaling algorithms for Kohn-Sham type models [25,38,24,6]. A stronger notion is the locality of the mechanical response, which is a fundamental premise underpinning the construction of interatomic potentials and of multi-scale algorithms such as hybrid QM/MM schemes [17] (here, it is termed "strong locality"). This latter category of locality is less well studied, the only result in this direction being the locality of non-selfconsistent tight binding models [14].
The aim of the present work is to establish the locality properties satisfied by the Thomas-Fermi-von Weizsäcker (TFW) model. Our main technical result to achieve this is the following pointwise stability estimate for the TFW equations, which establishes the locality of the electron response to changes in the nuclear configuration. Compared with [14] it is noteworthy that our result takes Coulomb interaction fully into account. Rigorous statements, under different conditions, are given in Theorems 3.4 and 3.5.
Theorem. For i = 1, 2 let m i ∈ L ∞ (R 3 ) represent nuclear charge distributions satisfying Let the corresponding ground state electron densities and electrostatic potentials, denoted by u i , φ i : R 3 → R, satisfy the TFW equations, Date: October 25, 2016. FQN is funded by the MASDOC doctoral training centre, EPSRC grant EP/H023364/1. CO was supported by ERC Starting Grant 335120 and the Leverhulme Trust through a Philip Leverhulme Prize. 1 Then there exists C, γ > 0 such that for all y ∈ R 3 (1. 1) In the remainder of the article we explore some of the consequences of this locality result: In Proposition 4.1 we obtain new estimates on finite-domain approximations which yield exponential decay of surface energies as well as an exponential convergence rate for the thermodynamic limit. In Corollary 4.2 we show that (1.1) gives rise to rigorous results that match, and substantially generalise, the Thomas-Fermi theory of impurity screening in metals [2,28]. In Theorem 4.3 we strengthen existing results on the neutrality of the TFW model [10]. In all these results, general (condensed) nuclear arrangements are treated.
A striking application of (1.1) is that it allows us to decompose energy into local contributions from which we obtain local site energy potentials: Given a countable collection of nuclei Y = (Y j ) j∈N ⊂ R 3 we construct an energy density E(Y ; x) which allows us to define the TFW energy Ω E(Y ; x) dx of an arbitrary volume Ω ⊂ R 3 in a meaningful way. This then motivates us to define site energies where (ϕ j ) j∈N is a smooth partition of unity of R 3 , which can be constructed in such a way that E j are permutation and isometry invariant and most crucially, E j are local in the sense that for some C, γ > 0. The rigorous statement of this result is given in Theorem 4.4. An analogous result has recently been proven for a tight binding model in [14]. This result not only gives a strong justification for the construction of classical shortranged interatomic potentials in metals, but in fact it allows us to treat the TFW mechanical response as if it emanated from such a classical potential. For example, (i) the analysis of the Cauchy-Born continuum limit [37] applies directly to the TFW model; and (ii) we can generalise in [13] the analysis of variational problems for the mechanical response to defects in an infinite crystal [19].
The remainder of this article is organised as follows: In Section 2 we recall the definition of the TFW model and summarise the relevant existing results. In Section 3 we state the main technical results, including the rigorous statement of the stability result (1.1). In Section 4 we present applications. Concluding remarks are made in Section 5, followed by the detailed proofs of the results in Section 6. Remark 1. We conclude the introduction with a remark about the analytical context of this work. The TFW equations for the electron density and potential is a coupled Schrödinger-Poisson system. Other systems of this class can be found in semiconductor physics [34,43,45]. We also note that Thomas-Combes estimates give conditions under which eigenfunctions of a Schrödinger operator decay exponentially [1,16]. While the results obtained for these systems are similar, the corresponding equations have different structure, hence the analytical techniques used to study them differ considerably.
The closest existing result to (1.1) we have found is [8,Theorem 4.6], which shows the exponential decay of the electron density away from the boundary of a crystal. Both (1.1) and [8,Theorem 4.6] utilise the uniqueness of the TFW equations to prove stability estimates. In Section 4.1, we use Proposition 4.1 to generalise [8,Theorem 4.6].
Acknowledgements. We thank Chen Huang for his valuable comments on impurity screening in the Thomas-Fermi model, and Eric Cances and Virginie Ehrlacher for helpful discussions about different analytical techniques for the TFW model.

(2.2)
Let m ∈ L 6/5 (R 3 ), m ≥ 0 denote the charge density of a finite nuclear cluster, then the corresponding TFW energy functional is defined, for v ∈ H 1 (R 3 ), by 3) The function v corresponds to the positive square root of the electron density. The first two terms of E TFW (v, m) model the kinetic energy of the electrons while the third term models the Coulomb energy. We remark that this definition of the Coulomb energy is only valid for smeared nuclei. We can rescale the energy to ensure C W = C TF = 1.
To construct the electronic ground state for an infinite arrangement of nuclei (e.g., crystals), we restrict admissible nuclear charge densities to m ∈ L 1 unif (R 3 ), m ≥ 0, satisfying (H1) sup The property (H1) guarantees that no clustering of infinitely many nuclei occurs at any point in space whereas (H2) ensures that there are no large regions that are devoid of nuclei. Let R n ↑ ∞ and define the truncated nuclear distribution m Rn = m · χ B Rn (0) , then the minimisation problem where C is independent of R n . Consequently, (2.5) implies that along a subsequence (u Rn , φ Rn ) converge to (u, φ). Passing to the limit in (2.4) yields the following result.
in the distributional sense. In addition, inf u > 0.

Remark 2.
A concise proof of uniqueness of the TFW equations is given in [7] under the assumption that m is smooth and hence u, φ ∈ W 1,∞ (R 3 ), which simplifies the earlier proof given in [12].
In the next section, we discuss results that can be obtained by generalising the proof of Theorem 2.1.

3.1.
Uniform regularity estimates. In the proof of our main results in the next section we employ regularity estimates that refine those of [12], and may be of independent interest.
Other than Proposition 3.1, our estimates rely on uniform variants of (H1)-(H2) and it turns out that (H2) may then be simplified without loss of generality; see Lemma 6.1 for more details. Given M, ω 0 , ω 1 > 0, let ω = (ω 0 , ω 1 ) and define the class of nuclear configurations As each nuclear distribution m ∈ M L 2 (M, ω) satisfies (H1)-(H2), Theorem 2.1 guarantees the existence of corresponding ground states (u, φ). We adapt the proof of existence of Theorem 2.1 to show that the uniformity in upper and lower bounds on m ∈ M L 2 (M, ω) yields regularity estimates and lower bounds on these ground states which are also uniform.
Proposition 3.1. For any nuclear distribution m : R 3 → R ≥0 , satisfying there exists (u, φ) solving (2.6) and satisfying u ≥ 0 and There exists c M,ω > 0 such that for all m ∈ M L 2 (M, ω) the corresponding ground state (u, φ) is unique and the electron density u satisfies Assuming higher regularity of the nuclear distributions implies higher regularity of the ground state. We therefore define, for k ∈ N 0 , Arguing by induction and applying the uniform lower bound (3.4) yields the following result.
3.2. Pointwise stability and locality. We now discuss our main result, a pointwise stability estimate for (2.6) which reveals a generic locality of the TFW interaction. To establish this result, we adapt the proof of uniqueness of the TFW equations in [12,7], specialising the class of test functions to for some γ > 0. Observe that e − γ|·| ∈ H γ for any 0 < γ ≤ γ.
We can generalise Theorem 3.4 to obtain higher-order pointwise estimates, but this requires both inf u 1 , inf u 2 > 0, hence we need to assume m 1 , m 2 ∈ M H k (M, ω) for some k ∈ N 0 . Theorem 3.5. Let k ∈ N 0 and m 1 , m 2 ∈ M H k (M, ω). Consider the corresponding ground states (u 1 , φ 1 ), (u 2 , φ 2 ) and define Then, there exist C = C(k, M, ω), γ = γ(M, ω) > 0 such that for any ξ ∈ H γ In particular, for any y ∈ R 3 , The existence and uniqueness of a corresponding ground state (u, φ) is guaranteed by [12,Theorem 6.10]. We believe that the arguments used to show [12,Lemma 5.5] and Theorem 3.4 can be adapted to show pointwise estimates similar to (3.8)-(3.9) when m 1 , m 2 satisfy (H1')-(H2') and that m 1 −m 2 is absolutely continuous with respect to the Lebesgue measure on R 3 , with a density belonging to L 2 unif (R 3 ). This result is not sufficient to consider the response of the ground state to a perturbation of point nuclei, though it may be possible to treat this using an approximation to the identity or by applying similar techniques. 4. Applications 4.1. Thermodynamic limit estimates. The following result provides an estimate for comparing the infinite ground state with its finite approximation, over compact sets, thus providing explicit rates of convergence for the thermodynamic limit. This is discussed in Remark 5.

4.2.
Locality of the Charge Response. The following result shows that the decay properties of the nuclear perturbation are inherited by the response of the ground state. Consider the corresponding ground states (u 1 , φ 1 ), (u 2 , φ 2 ) and define Remark 6. For some of our comparison results, we require only m 1 ∈ M L 2 (M, ω) but impose weaker assumptions on m 2 . This would not generalise Corollary 4.2 since any of the decay assumptions in (1-3) already imply that m 2 ∈ M L 2 (M, ω ′ ) for some ω ′ .
Remark 7. The estimate (4.3) can be used to study the full non-linear response of the ground state to a nuclear impurity. We compare this to the results from the Thomas-Fermi (TF) [2,28,39] and TFW [18,35,40,29] theories of screening. Consider a nuclear arrangement m 1 ∈ M L 2 (M, ω) and model a nuclear impurity at the origin with positive charge Z by Zη(x), where η ∈ C ∞ c (R 3 ), η ≥ 0 and η = 1. Then define the perturbed system by m 2 = m 1 + Zη ∈ M L 2 (M 1 , ω 1 ) and consider the corresponding TFW ground states (u 1 , φ 1 ) and (u 2 , φ 2 ), respectively. From (4.3) of Corollary 4.2 it follows that We now compare (4.6) with existing results from the TF and TFW theories of screening. These models consider the formal linear response (n, V ) of the electron density and potential to a nuclear impurity at the origin, modelled by the Dirac distribution Zδ 0 , in a uniform electron gas. In both models, V satisfies the linear equation where k s > 0 is a material-dependent constant called the inverse screening length. In the TFW model, V and n satisfy [18,40,29,35] V where α ∈ R, β ∈ C satisfy 0 < |β| < α. The constants α, β depend on the material and the coefficient C W , which appears in the definition of the TFW energy (2.3). There is a critical value of C W below which β > 0 and above which β is complex, the latter case introduces oscillations in the potential and electron density. In either case, both the TF and TFW models exhibit screening due to the presence of the exponential term appearing in (4.7)-(4.8).
The lack of a factor of the form 1 |x| in (4.6) can be attributed to using a smeared nuclear description for the impurity as opposed to a point description in (4.7)-(4.8). Other than this, the similarity of (4.6) to (4.7) suggests that the constant γ in (4.6) may be interpreted as the inverse screening length. In this paper we show there exists γ > 0 satisfying (4.6), however we do not provide any estimates for its value.
The estimate (4.6) shows that screening occurs in the TFW model, without any approximations made to the model and without any restrictions on the nuclear configurations (other than (H1)-(H2)). It should be noted that although (4.6) agrees with existing results from the TF theory of screening, in metals often the effects of screening are weaker. For metals, instead of an exponentially decaying screening factor, Friedel oscillations are observed [22,33,27]. In this case, the screening factor behaves as |x| −r f (|x|), where f : R ≥0 → [−1, 1] is an oscillating function and the decay rate r > 0 depends on the Fermi surface of the metal. The generic exponential screening factor in (4.6) demonstrates that the TFW model significantly overscreens charges. 4.3. Neutrality of defects. An immediate consequence of Corollary 4.2 is the neutrality of nuclear perturbations in the TFW equations. This result applies to all nuclear configurations belonging to M L 2 (M, ω). In particular Theorem 4.3(3) strengthens the result of [10], which requires m 1 − m 2 ∈ L 1 (R 3 ) ∩ L 2 (R 3 ) and thus excludes typical point defects; see Remark 8 for more details.
Remark 8. In a forthcoming article [13], we construct a variational problem to study the response of a crystal due to a local defect, using the TFW energy. Arguing as in [19], we shall show that any minimising displacement decays away from the defect at the rate |x| −2 , which corresponds to case (2) with r = 2. In this case (4.10) only provides a uniform bound for the charge as opposed to a decay estimate. However, as r > 3/2 the global neutrality result (4.11) holds for the relaxed system. The neutrality estimates of Theorem 4.3 strengthen those of [10] in the following ways. Firstly, our result considers a perturbation of a general nuclear arrangement as opposed to a perfect crystal. This allows us, in [13], to consider the response of extended defects such as dislocations. In addition, we only require that the nuclear perturbation belongs to L 2 (R 3 ), which we prove rigorously in [13], whereas in [10] the nuclear defect is assumed to lie in L 1 (R 3 ) ∩ L 2 (R 3 ), which fails for typical point defects.

Energy locality.
We now show that the locality result, Theorem 3.5, can be used to describe the energy contribution of each individual nucleus. In effect, we will derive a site energy potential for the TFW model, which has the surprising consequence that, for the study of mechanical response, TFW can be treated as a classical short-ranged interatomic potential. Our result gives credence to the construction of interatomic potentials and the assumption of strong locality used in hybrid quantum mechanics/molecular mechanics (QM/MM) simulations [17].
Let η ∈ C ∞ c (B R 0 (0)) be radially symmetric and satisfy η ≥ 0 and R 3 η = 1 describe the charge density of a single (smeared) nucleus, for some fixed R 0 > 0. For any countable collection of nuclear coordinates Y = (Y j ) j∈N ∈ (R 3 ) N , let the corresponding nuclear configuration be defined by (4.12) A natural space of nuclear coordinates, related to the M L 2 space is This space contains many condensed phases, such as crystals containing point defects, dislocations and grain boundaries. It does not include arrangements with arbitrarily large voids such as surfaces or cracks. However, as the TFW model for surfaces has been established [7], it may be possible to obtain locality estimates for surfaces and cracks using the TFW model. We remark that there exists For any Y ∈ Y L 2 (M, ω) there exists a unique ground state (u, φ) corresponding to m = m Y . Naively, we might define the energy stored in a region Ω ⊂ R 3 by however, difficulties arise due to the fact that (m − u 2 ) * 1 |·| is not well-defined. Instead, we give two alternative definitions for the energy density for an infinite system: . Suppose now that Ω ⊂ R 3 is a charge-neutral volume [44], that is, if n is the unit normal to ∂Ω, then ∇φ · n = 0 on ∂Ω. Recalling from (2.6b) that In particular, for finite neutral systems m Rn = m · χ B Rn (0) , where m ∈ M L 2 (M, ω) and (u Rn , φ Rn ) denoting the corresponding solution, the following energies agree on Ω = R 3 We prove this claim is in Remark 12 in Section 6. Thus, we have derived two energy densities, E 1 , E 2 , which are meaningful and well-defined also for infinite configurations. In order to define site energies, we require a partition of R 3 . For each j ∈ N let ϕ j (Y ; ·) ∈ C 1 (R 3 ), ϕ j (Y ; ·) ≥ 0 satisfying the following conditions: there exist C, We propose a canonical construction of such a partition in Remark 9 below. Given a family of partition functions satisfying (4.19), we can define site energies for i = 1, 2. A consequence of Theorems 3.4 and 3.5 is that E i j (Y ) are local: their dependence on the environment of nuclei decays exponentially fast. This is made precise in the following theorem.
The derivative ∂ Y k E i j can be interpreted as the contribution of the atom at Y k to the force on the nucleus at Y j . In addition, we show in Section 6.4 that these site energies generate the correct total force Remark 9. Two further canonical requirements on a site energy potential are permutation and isometry (rotation and translation) invariance. This can be obtained as follows: If the partition (ϕ j ) j∈N is permutation invariant, that is, for any bijection P : then so are the site energies, then the site energies are also isometry invariant, . Both statements are proven in Lemma 6.9.
A canonical class of partitions satisfying (4.19) as well as (4.23), (4.24) can be constructed as follows: Let ϕ ∈ C 1 (R 3 ), ϕ ≥ 0, be radially symmetric and satisfy For example, this holds for ϕ(x) = e − γ|x| 2 and for standard mollifiers with sufficiently wide support. Then, for j ∈ N, we can define .

(4.27)
It is easy to see that this class of functions are well-defined and satisfies all requirements.
Remark 10. Alternative constructions of energy partitions include Bader volumes and charge-neutral volumes [4,44,32]. Bader volumes partition space into regions such that the flux of the electron density on the boundary is zero, while charge-neutral volumes are defined so that each region has zero charge. The construction of these volumes is not unique, like our definition of a partition. Bader volumes were developed as a means to define atoms within molecules [4].
With this in mind, using a partition we may assign a portion of the electron density to each nucleus in the system. We refer to a nucleus paired with it's associated partition of the electron density as an effective atom. Due to the screening that occurs in the TFW model, the interaction of two effective atoms decays exponentially as the distance between the nuclei grows. In comparison, the interaction of two neutrals atoms separated by a sufficiently large distance r in the TF model has been shown to decay at the rate r −7 [9]. This suggests that due to the overscreening of the TFW model, the interaction of the effective atoms is considerably weaker than is realistic. However, for the purpose of simulating quantum systems, in particular applying the strong locality principle [17], the weak long-range interaction of the TFW model is a desirable property. Remark 11. The estimate shown in Theorem 4.4 is a theoretical result which can be applied to simulate defective crystals, though we do not pursue this. Locality estimates have been established for the tight-binding model and subsequently used to construct QM/MM hybrid methods [14,15].

Conclusion and Outlook
The two main results of this work, Theorems 3.4 and 3.5, are stability and exponential locality estimates for the TFW model, which apply to general condensed nuclear configurations.
We have demonstrated in Section 4 that it can be used to extend and strengthen a range of existing results on the TFW model. A particular strength of our results is that they apply to general nuclear configuration in M L 2 (M, ω), whereas the previous analyses of the TFW model have focused on (near-)crystalline arrangements or the homogeneous electron gas. This generality will be valuable when exploring the consequences of our analysis for studying models for the mechanical response problem in [13], where we generalise [19] to electronic structure models.
A further application, that we will develop in a forthcoming work is a study of the Yukawa potential as a model approximation [36]. Adapting Theorems 3.4 and 3.5 we can consider the difference between the Coulomb and Yukawa ground states for a given nuclear configuration and prove uniform error estimates in terms of the screening parameter in the Yukawa model.
Two key difficulties in the analysis of electronic structure models are (i) the exchange and correlation of electrons due to the antisymmetry of the electronic wavefunction; and (ii) the interaction of charged particles (positive nuclei and negative electrons) via the longrange Coulomb potential. Since the TFW model is orbital-free it does not account for (i), however it fully incorporates Coulomb interaction. In this regard it is perhaps surprising that the TFW model satisfies the extremely generic locality property we obtained in Theorems 3.4 and 3.5.
The Hartree-Fock and Kohn-Sham models take both effects into account and whether these models permit a similarly strong notion of locality is an open problem. It is clear, however, that such results cannot be obtained in the generality that we considered in the present paper. Since charged defects exist in the reduced Hartree-Fock model [11] and as locality implies neutrality, this suggests that a locality property cannot hold for general condensed phase arrangements in the reduced Hartree-Fock model, which is the simplest model in the Hartree-Fock/Kohn-Sham class.

Proofs
This section contains the proofs of the main results. Proofs of results in Sections 3.1, 3.2 and 4 are found in Sections 6.1, 6.2 and 6.3 respectively.
The following is a preliminary result used in the construction of the space M L 2 (M, ω).
Proof of Lemma 6.1. Clearly, (6.1) implies (H2), so suppose m satisfies (H2), then there exists Now define ω 1 := ω 0 R 3 2 ≥ 0, then it follows from (6.2) that (6.1) holds for all R > 0. 6.1. Proofs of Uniform Regularity Estimates. The following lemma features in the proofs of both the existence and uniqueness of the TFW equations and is found in [12].
, then define the elliptic operator L = −∆ + a. Suppose that there exists u ∈ H 1 loc (R 3 ) satisfying u > 0 and Lu = 0 in distribution. Then, the operator L is non-negative, that is for all The proof is shown in [12] but is included here for completeness.
Proof of Lemma 6.2. Let R > 0 and define Ω = B R (0) and consider L as an operator on L 2 (Ω) with domain H 2 (Ω) ∩ H 1 0 (Ω). Then L is a self-adjoint operator with compact resolvent hence has a purely discrete spectrum. Since a ∈ H 1 loc (R 3 ) it follows that the smallest eigenvalue λ 1 (Ω) is simple and has a positive eigenfunction v Ω ∈ H 1 0 (Ω) [23,Theorem 8.38]. In addition, by standard elliptic regularity v Ω ∈ H 3 (Ω) ֒→ C 1,1/2 (Ω) [20] and solves Testing this equation with u and using integration by parts, we obtain As v Ω > 0 on Ω and v Ω vanishes over ∂Ω, it follows that ∂v Ω ∂n ≤ 0. It follows that the left-hand side of (6.4) is non-negative, hence λ 1 (Ω) ≥ 0. As this holds for Ω = B R (0), for any R > 0, we deduce that for all ϕ ∈ C 1 c (R 3 ) ϕ, Lϕ ≥ 0. Using that a ∈ L ∞ (R 3 ) and the density of , it follows that for all ϕ ∈ H 1 (R 3 ) ϕ, Lϕ ≥ 0. We now show uniform estimates for finite systems corresponding to truncated nuclear distributions. This result is essentially [12,Proposition 3.5], however as we require uniform regularity estimates, we provide a complete proof.
and R n ↑ ∞, then define the truncated nuclear distribution m Rn = m·χ B Rn (0) . The unique solution to the minimisation problem Rn ). which satisfy the following estimates, with constant C independent of R n : Proof of Proposition 6.3. If m ≡ 0, then for all R n , clearly u Rn = φ Rn = m Rn = 0 satisfies (2.4) and (6.7)-(6.8).
If m ≡ 0, then there exists a constant R 0 ≥ 0 such that R n ≥ R 0 ensures that For each R n , consider the minimisation problem The constraint R 3 v 2 = R 3 m Rn ensures the system is charge neutral, and by [31,Theorem 7.19] there exists a unique non-negative minimiser Here θ Rn > 0 is the Lagrange multiplier associated with the charge constraint (6.10) [31,12]. Define φ Rn : so we can express (6.9) as the Schrödinger-Poisson system (2.4) .
Using that u Rn , φ Rn + θ Rn are continuous and decay at infinity, by arguing as in [41], there exists a universal constant C S > 0, independent of the nuclear distribution, satisfying 0 < θ Rn ≤ C S , (6.14) 10 9 u 4/3 As u Rn ≥ 0, from the Solovej estimate (6.15) we obtain the uniform lower bound We aim to show a uniform upper bound for φ Rn , which together with (6.15) will yield the uniform estimate which is independent of R n . If φ Rn is non-positive, then (6.17) holds as Rn is non-zero at some point in R 3 . By (6.11) φ Rn is a continuous function that converges to a negative limit at infinity, Without loss of generality, we assume x Rn = 0. We now show that u Rn > 0 on R 3 , arguing by contradiction. Suppose that there exists z ∈ R 3 such that u Rn (z) = 0. Since u Rn is a non-negative, continuous function decaying at infinity, there exists y n ∈ R 3 such that u Rn (y n ) = sup Let R > |y n − z|, then by the Harnack inequality [42], we infer u Rn (x) = u Rn (z) = 0, so u Rn ≡ 0. This contradicts the charge constraint (6.10) Rn − φ Rn is a non-negative operator. Choose ϕ ∈ C ∞ c (B 1 (0)) satisfying 0 ≤ ϕ ≤ 1, ϕ = 1 on B 1/2 (0), R 3 ϕ 2 = 1 and R 3 |∇ϕ| 2 =: c ϕ , then for y ∈ R 3 , define ϕ y ∈ C ∞ c (B 1 (y)) by ϕ y = ϕ(· − y). As L Rn is non-negative (6.3) implies which can be re-arranged and expressed using convolutions as 5 3 u 4/3 We estimate the first term using (6.5) For the second term, using the convexity of t → t 3/2 for t ≥ 0 and the fact that ϕ 2 = 1, applying Jensen's inequality and (6.19) we deduce Rn (x − y)ϕ 2 (y) dy Combining the estimates (6.20) -(6.22) we conclude that −∆ φ Rn * ϕ 2 + φ Rn * ϕ 2 − c ϕ 3/2 + ≤ C 0 M. By (6.11), as φ Rn is a continuous function that converges to a negative limit at infinity, φ Rn * ϕ 2 also shares these properties. Now consider the set it follows that S is open and bounded and that φ Rn * ϕ 2 − c ϕ = 0 on ∂S. Observe that the constant function h = (C 0 M) 2/3 satisfies where , and by estimating (2.4b) directly, that Rn . As 0 ≤ ϕ ≤ 1 and ϕ = 1 on B 1/2 (0), then Using a change of variables, (6.24) can be expressed as and suppose that for all t ∈ (1/4, 1/2) which gives a contradiction, hence for some t ∈ (1/4, 1/2) We now construct an upper bound for φ + Rn as follows. Let φ 1 satisfy −∆φ 1 = 0 in B t (0), As φ 1 is harmonic, it satisfies the mean value property Then consider the Dirichlet problem By Lax-Milgram, this has a unique weak solution φ 2 ∈ H 1 0 (B t (0)). By standard elliptic regularity theory [20] φ 2 ∈ H 2 (B t (0)) ֒→ C 0,1/2 (B t (0)) and hence by the maximum principle φ + Rn ≤ φ 1 + φ 2 , in particular (6.26)-(6.27) imply where the right-hand side is independent of R n . Combining this with the lower bound (6.16) and the Solovej estimate (6.15), we obtain the estimate (6.17) independently of both x, r and R n . Using (6.17) and (6.28), we now obtain uniform local estimates for the right-hand side of (6.13) Rn + φ Rn u Rn by Consequently, for any x ∈ R 3 , the elliptic regularity estimate [20] gives Applying a similar argument to estimate the right-hand side of (2.4b) −∆φ Rn = 4π(m Rn − u 2 Rn ) yields (6.8) ). Using that φ Rn ∈ H 2 unif (R 3 ) and arguing as in (6.30), we obtain the desired estimate (6.7) u Rn H 4 unif (R 3 ) ≤ C(1 + M 15/4 ). We remark that while the constants appearing in the final estimates (6.7)-(6.8) depend on c ϕ , they are independent of M.

Remark 12.
We now justify the claim that for finite and neutral systems and for Ω = R 3 , the three energies shown in (4.18) agree. Recall (6.12), which shows that the Coulomb energy can be expressed as so it follows that the energies defined in (4.18) agree for Ω = R 3 .
We now discuss passing to the limit in (2.4) to obtain regularity for the infinite system.
such that along a subsequence u Rn , φ Rn converges to u, φ, weakly in H 2 (B R (0)), strongly in H 1 (B R (0)) for all R > 0 and pointwise almost everywhere. It follows from the pointwise convergence that u ≥ 0 and Passing to the limit of the equations (2.4) in distribution, we find that the limit (u, φ) solves Arguing as in (6.7)-(6.8), we deduce that the desired estimates (3.2)-(3.3) hold we show that this contradicts the assumption that for all m ∈ M L 2 (M, ω) and R > 0 It follows from (6.31) that there exists m n ∈ M L 2 (M, ω) with corresponding solution (u n , φ n ) and x n ∈ R 3 such that for all n ∈ N u n (x n ) ≤ 1 n .
Recall the uniform estimates (6.7)-(6.8) from Proposition 3.1 It follows that the sequence of functions u n (·+x n ) converges uniformly to zero on compact sets. Consider the ground state (u n , φ n ) corresponding to the nuclear distribution m n .
Recall that φ n solves the following equation in distribution −∆φ n = 4π m n − u 2 n . (6.32) We translate the system and then pass to the limit in (6.32) as n tends to infinity. To do this, we use the following estimates, which are translation invariant: . It follows that, up to a subsequence, φ n (· + x n ) converges to φ, weakly in H 2 (B R (0)), strongly in H 1 (B R (0)) for all R > 0 and pointwise almost everywhere. Moreover, m n (· + x n ) converges to m, weakly in L 2 (B R (0)) for all R > 0. By applying the Lebesgue-Besicovitch Differentiation Theorem [21] we deduce that m ∈ M L 2 (M, ω). Passing to the limit in Arguing as in [12, Theorem 6.10], we show that for all R > 0 As m ∈ M L 2 (M, ω), this leads to the contradiction that for all R > 0 To show (6.34) choose ϕ ∈ C ∞ c (B 2 (0)) such that 0 ≤ ϕ ≤ 1 and ϕ = 1 on B 1 (0). Let R > 0, then testing (6.33) with ϕ(·/R) gives The left-hand side can be estimated by where the constant C > 0 is independent of R. As m ≥ 0, from (6.35) we obtain (6.34) . We now show the induction step. Suppose the result is true for k ∈ N 0 , then consider m ∈ M H k+1 (M, ω) ⊂ M H k (M, ω), so by the induction hypothesis the corresponding solution (u, φ) satisfies We remark that as 0 < c M,ω ≤ u ≤ C(M) and u ∈ H k+4 unif (R 3 ), it follows that for all r ∈ R, by standard elliptic regularity theory [20] for any x ∈ R 3 We use an identical argument together and apply the estimate (6.38) to deduce Combining (6.38) and (6.39) we obtain the desired estimate , ω , which completes the proof of (6.36) by induction.
6.2. Proofs of Pointwise Stability Estimates. To prove Theorems 3.4 and 3.5, we adapt the proof of uniqueness of the TFW equations, shown in [12,7]. Due to the length of the argument, we shall prove several intermediate results. Before showing these results, we outline the structure of the proof.
Throughout the remainder of the paper we use the notation By treating the coupled system of equations as a linear system and by exploiting the coupling between the electron density and electrostatic potential arising from the Coulomb energy term of the TFW functional, we obtain the following initial estimates To control the ψ-dependence on the right-hand side of (6.41), we require an estimate of the form To remove the term (w 2 + ψ 2 )|∇ξ| 2 on the right-hand side, we simply restrict from ξ ∈ H 1 to a narrower class of test functions, In order to show (6.42), we adapt the argument used in [7] . At the same time, since the equations for (w, ψ) hold pointwise, we obtain additional estimates for ∆w, ∆ψ.
Clearly Lemmas 6.4 and 6.5 hold also under the assumption (B) since (B) implies (A), with M ′ = C(M). In the case (B) where m 1 , m 2 are both uniformly bounded below and have higher regularity, arguing as in Corollary 3.3 and Lemma 6.5, we obtain improved estimates for w and ψ.
Observe that in Case (B), M ′ = C(M). Due to this, we omit the dependence of M ′ in the constants that appear in the following lemmas, whenever we assume (B) holds. Lemma 6.6. Suppose that either (A) or (B) holds, then there exist In particular, for any y ∈ R 3 , We remark that in the following proofs, all integrals are taken over R 3 , unless stated otherwise.
Combining (6.45) and (6.79), we obtain the desired estimate (6.46) . The argument used to show (6.70) holds for k ′ ≤ 2, so for ξ ∈ H 1 Then, as (6.75) holds with j 1 ≤ 2, applying this and (6.70) for k ′ = 0 yields Similarly, the argument used to show (6.71) holds for k ′ = 0, to give Finally, combining (6.80)-(6.81) and applying (6.44) from Lemma 6.5, we obtain the desired estimate (6.45) with k = 0 The argument used in Case 1 holds for k = 0 to show the desired estimate (6.46) We have now established all technical prerequisites to prove Theorems 3.4 and 3.5.
Next, we now prove Corollary 4.2 as a direct consequence of Theorems 3.4 and 3.5.
Under the assumptions of Theorem 3.4 with k = 0, other than applying Theorem 3.4 instead of Theorem 3.5, the proof is identical.
We turn to the proofs of the charge-neutrality estimates.
6.4. Proof of energy locality. To prove Theorem 4.4, we first establish the existence, uniqueness and regularity of the solutions to the linearised TFW equations.
We now verify that the limiting functions are independent of the choice of sequence. First, observe that by passing to the limit as h n → 0 in the equations it follows that (u, φ) solve the linearised TFW equations (6.96) pointwise, −∆u + 35 9 u 4/3 − φ u − uφ = 0, Clearly m is independent of the sequence h n . Applying [7,Corollary 2.3], it follows that the (u, φ) is the unique solution to the linearised system (6.96), hence is independent of the sequence (h n ). It then follows that u h −u h , φ h −φ h converge to u, φ as h → 0 as stated above.
We are now in a position to prove Theorem 4.4.
Proof of Theorem 4.4. We will repeatedly use the fact that there exists C, γ > 0 such that, for all h ∈ [0, h 0 ], p ∈ [1, 2], which is a consequence of the uniform bounds on m h , φ h and of (6.86). Further, we require that there exist C, γ > 0 such that, for j ∈ N, h ∈ (0, h 0 ], x ∈ R 3 , which follows directly from (4.19c). For i = 1, 2 and j ∈ N, consider the difference We wish to show that the limit of (6.113) exists as h → 0 to obtain To show (6.115), consider the difference It follows from (6.118) and pointwise convergence of u h , ∇u h , φ h to u, ∇u, φ and to u, ∇u, φ, m as h → 0, that (6.115) holds Applying (6.94) to (6.118) yields Combining (6.119) and (4.19b), we deduce hence by (6.111) and the Dominated Convergence Theorem, It follows from (6.120) and (6.111) that (6.122) It remains to show that (6.113) converges using (6.111) and (6.112). As ϕ j (Y ; x) is differentiable with respect to Y k , for all x ∈ R 3 and combining (6.117) with (6.112) implies hence by (6.111) and the Dominated Convergence Theorem, Combining (6.122) and (6.123) yields the desired estimate (4.21). The second case, using E 2 instead of E 1 , is analogous.