Global Existence for the N Body Euler-Poisson System

In this paper we investigate the problem of multiple expanding Newtonian stars that interact via their gravitational effect on each other. It is clear physically that if two stars at rest are separated initially, and start expanding as well as moving according to the laws of Newtonian gravity, they may eventually collide. Thus, one can ask whether each star can be given an initial position and velocity such that they can keep expanding without touching. We show that even with gravitational interaction between the bodies, a large class of initial positions and velocities give global-in-time solutions to the N Body Euler-Poisson system. To do this we use a scaling mechanism present in the compressible Euler system shown in [63] and a careful analysis of how the gravitational interaction between stars affects their dynamics.


Introduction
The compressible Euler-Poisson system for an inviscid, isentropic, ideal gas, acting under the influence of its own gravity, in its most basic form, is given by where ρ, u, p, and φ are the fluid density, velocity, pressure, and gravitational potential respectively. Newton's gravitational constant G will be set to 1 here and hereafter in this work.
Remark 1.1. Investigating the properties of stars modelled as bodies of fluids is a classical problem in astrophysics, with a long history. In the case of one body of incompressible fluid moving under its own gravity, equilibrium states have been studied by the likes of Newton, Maclaurin, Jacobi, Poincaré, and others. See, for example, [8] for a summary of these results, and further context regarding this problem.
Remark 1.2. Note that if instead we had ∆φ = −4πρ, φ would be referred to as the electrostatic potential, and the above would be a toy model of a plasma instead of a Newtonian star. See [26,27] for more details on the mathematical theory of plasmas.
In addition to the basic system (1)-(3), we will always work with a polytropic equation of state between the pressure p and the density ρ: for some γ > 1, the adiabatic index. The system (1)-(4) is one of the most fundamental and commonly used models of a Newtonian star ( [7,79,4]). In addition, the equation of state guarantees the system is not underdetermined. A famous class of solutions to (1)-(4) are the so called Lane-Emden stars, obtained by looking for radially symmetric, time independent density profiles, and vanishing velocity: (ρ(t, x), u(t, x)) = (ρ(r), 0).
In the range 6/5 < γ < 2, it is known that there always exists a solution with finite mass and compact support, whilst for γ = 6/5 there exists a solution with finite mass and infinite support ( [7,79,4,35]). The fundamental problem
Remark 1.3. Here κ is introduced as an index variable ranging from 1 to N . Throughout this work, indexing the star bodies will be the one and only use of κ, and conversely, indexing of the separate stars will only be done with κ.
In particular there will be no summation convention over κ.
Analogously to (1)-(3), the ρ κ , u κ , and p κ are the densities, velocities, and pressures for each separate star. As mentioned in (4), we have a polytropic equation of state for each κ given by where γ > 1 is constant, and the same for each κ.
Remark 1.4. It is not crucial for each γ to be the same for each equation of state, i.e. γ could depend on κ. However, it simplifies the analysis at certain points.
Comparing (1)-(4) to (6)- (12), we see that in the latter system, for each fixed t, the continuity and momentum equations for each κ are only required to hold on the support of ρ κ (t, ·), namely Ω κ (t), instead of on all of R 3 as in the former system. In addition, the support, and therefore its boundary is now a dynamic object whose evolution we track in our framework, making this model a vacuum free boundary problem.
The ρ κ are non-negative functions ρ κ : [0, ∞) × R 3 → R ≥0 , with supports supp ρ κ (t, ·) = Ω κ (t). The boundary of Ω κ (t) is given by ∂Ω κ (t). We also have the initial domains and their boundaries Ω κ := Ω κ (0), (13) Since we are looking at stars that are separated from each other initially, physically we must have that the Ω κ are mutually disjoint and compact. The cumulative density function ρ : [0, ∞) × R 3 → R ≥0 is defined by and φ is the cumulative gravitational potential. For each fixed t, the domain of the velocity field at time t is given by the support of the density at time t: x → u κ (t, x).
Finally, the quantity V(∂Ω κ (t)) is the outward normal velocity of ∂Ω κ (t) and n κ is the outward unit normal on ∂Ω κ (t).
To construct a robust well-posedness theory for the moving vacuum boundary Euler-Poisson system, it is crucial to include the physical vacuum boundary condition. For each κ, the speed of sound c κ is given by Then the physical vacuum boundary condition reads −∞ < ∂c 2 κ ∂n κ ∂Ωκ < 0, (17) wheren κ is the outward unit normal on ∂Ω κ . Condition (17) implies that for some constant C which in turn implies thatc κ is only Hölder continuous of exponent 1/2 at the boundary Ω κ . We note that as well as being a key component for the local theory of our system, the physical vacuum condition also naturally occurs in Lane-Emden stars [35]. Liu [48] gave an argument suggesting that the physical vacuum condition was natural in the context of the Euler system with damping. In general, the physical vacuum condition, and the role it plays in understanding the interaction of fluids with vacuum regions has been a subject of great interest, see [49,50,51,11,12,13,36,37,38,39,54,69,30,31]. When we view the problem on R 3 , any sufficiently regular nontrivial spherically symmetric solution with compact initial data has a finite time of existence [57,56], with the same holding for any nontrivial solution for the Euler system without gravitation [58]. Results on singularity formulation for the Euler system have also been established by Sideris [70]. The solutions in [57,56] cannot satisfy the physical vacuum boundary condition (17), as the initial speed of soundc is continuously differentiable on the boundary.
In the case of the physical vacuum free boundary, well-posedness theories for the compressible Euler system were developed by Coutand and Skholler [13], and Jang and Masmoudi [39] independently. In our work, we adapt the techniques of Jang and Masmoudi. This theory is readily adapted to the Euler-Poisson system for the case N = 1 as the potential term ρ∇φ is lower order, with respect to derivative count, to the top order pressure term ∇p. We shall see in Appendix A that the local-in-time theory for N ≥ 2 does not require much adjustment in the theory either.
It is natural to look for expanding solutions to the Euler-Poisson system. Indeed, it has been shown that globalin-time solutions to the Euler system must have supports whose diameters grow at least linearly in time, with the same being true for solutions to the Euler-Poisson system with positive energy when γ ≥ 4/3, see [28,71]. Note that equilbrium solutions like the ones mentioned in the introduction are not considered in the class of time dependent solutions.
In the vacuum free boundary setting for the Euler system without gravitation, a finite parameter family of expanding global-in-time solutions has been found by Sideris [71,72], relying on an affine ansatz on the Lagrangian flow map, allowing him to reduce the problem to one of solving a system of ODEs. The corresponding solutions are a finite dimensional family of compactly supported, expanding stars. Using an affine ansatz to construct solutions for compressible fluid flow is a technique that goes back to Ovsiannikov [62] and Dyson [16]. Hadžić and Jang [30] showed nonlinear stability of the solutions constructed in [72] for the range γ ∈ (1, 5/3], after which Sideris and Shkoller [67] established the corresponding result for γ > 5/3. In the nonisentropic case, nonlinear stability of the Sideris solutions was established by Rickard, Hadžić, and Jang [65]. In addition, Hadžić and Jang [31] also showed that small perturbations of the Sideris affine solutions give rise to solutions to the vacuum free boundary Euler-Poission system, utilising a scaling structure in the system that allowed them to construct these solutions under the assumption of small densities. In general their methods require γ ∈ (1, 5/3). In this range, the gravitational potential terms are subcritical with respect to the pressure term, meaning that the pressure term will dominate the dynamics of the star. This is part of what allowed them to peturb solutions of the Euler system to find their solutions of the Euler-Poisson system. More specifically, their range of γ was restricted to {1 + 1 n |n ∈ Z ≥2 } ∪ (1, 14/13). Note that these restrcitions, apart from being a subset of (1, 5/3), are largely technical and related to estimating the gravitational potential terms.
The author, along with Hadžić and Jang [63], then showed that one can construct a family of global-in-time expanding solutions to the vacuum free boundary Euler and Euler-Poisson systems without appealing to an affine type ODE reduction. The main tool was utilising a scaling structure in the nonlinearity of the momentum equation that produced a stabilising effect. This, as well as the scaling exhibited in [31], was used to construct solutions with small densities. Analogously to [31], the range of γ for the Euler-Poisson system is {1 + 1 n |n ∈ Z ≥2 } ∪ (1, 14/13). Despite the smallness condition, the initial densities in [63] have a wide class of possible profiles, whereas those associated with the solutions exhibited in [72] have a very specific form.
In the absence of a free boundary, global-in-time expanding solutions have been found by Serre [68], Grassin [24], and Rozanova [66]. Serre found solutions by perturbing around linear velocity profiles, an idea which Grassin generalised. Rozanova showed the existence of related solutions without the small density assumption made by Grassin. With the polytropic equations of state in (12), our system becomes, for κ = 1, 2, . . . , N , We shall refer to the system (19)- (24) with the physical vacuum condition (17) as EP(N, γ). If we set in the vicinity of the initial vacuum boundary ∂Ω κ . The quantity N κ=1 w κ is proportional to the enthalpy of the system and the {w κ |κ = 1, 2, . . . , N } will play a very important role in our analysis. For future simplicity, define so thatρ κ = w α κ . By contrast to the classical N Body problem, the stars described by our system are subjected to tidal forces which deform the geometry of their supports, and hence they can not be idealised as point particles. The classical problem consists of looking for solutions to the system of N point particles interacting with each other via Newtonian gravity: where the x i (t) are the particle positions at time t, the m i are the particle masses,x i are their initial positions, and v i are their initial velocities. Much like the subject of stellar dynamics, the N body problem has a long and rich history; for more on the history and context, see [73,32] and references within. In our context, as mentioned previously, McCann [59] exhibited solutions to the stationary compressible Euler-Poisson system that corresponded to a rotating binary star sysem. Miao and Shahshahani [60] showed that in the case of the 2 body problem for the incompressible free boundary Euler-Poisson system, one can start with initial configurations that would correspond to hyperbolic orbit in the point particle case, and obtain bounded orbits. The reason is that the fluid bodies will naturally deform due to their gravitational effects on each other, and this results in a loss of energy.
In the classical N Body problem (28), one can ask what happens to the dynamics of each particle if we start them far away from each other (in some suitable sense), and also point their initial velocities away from each other. We might expect that if the initial distances are large enough, with initial velocities pointing in suitable directions, their dynamics will essentially decouple, and each will travel on a path of constant velocity up to small error. We will show in our work that this type of initial configuration can lead to global-in-time solutions even in the case of the N Body of Euler-Poisson system.
Explicitly, we exhibit an open set (in a suitable topology) of initial positions and velocities that lead to globalin-time solutions of EP(N, γ), with each star asymptotically behaving like an expanding star moving with constant velocity. This is further discussed in Sections 3.1, 6.1. Here we state a rough version of our theorem. Remark 1.6. The restriction of possible values of γ in our work is for the same reasons as in [31,63], discussed previously in this introduction.
The plan of the paper is as follows. In Section 2 we formulate our problem in Lagrangian variables, a technique commonly used to study the physical vacuum free boundary compressible Euler(-Poisson) system. In Section 3, we discuss the properties of scaling in the Euler-Poisson system, and how we use these properties to find global-in-time solutions. In Section 4, we fix general notation, and in Section 5 we define our energy spaces and higher order energy and curl functions. In Section 6 we state our main result precisely, and list our a priori assumptions. Sections 7, 8, and 9 are devoted to estimates of the graviational potential, curl estimates, and energy estimates respectively. Finally, in Section 10, we prove our main result, Theorem 6.5.

Lagrangian Formulation
For simplicity, we will assume that the initial domains Ω κ are closed balls of radius 1 with centresx κ . We introduce a reference domain: where B 1 is the closed unit ball centred at 0 in R 3 . For κ = 1, 2, . . . , N , it is clear that Ω κ = Ω +x κ .
Remark 2.1. The techniques we use to construct global-in-time solutions in this paper work equally well when each Ω κ is a diffeomorphism of a translate of the closed unit ball, close enough to the identity in some norm, with each translation such that the distance bwteen every pair of domains is sufficiently large.

Flow Map
In this section we introduce our ansatz for η κ . However, we first give a motivation for the form of our ansatz. In [63] the author alongside Hadžić and Jang found that in the case of N = 1, solutions of the form for the Euler and Euler-Poisson systems were global-in-time solutions for small enough θ (with size measured in a suitable function space), as well as small enough initial density. The suggested form of η in (48) is exactly what allows us to take advantage of the nonlinear scaling structure in the Euler, and Euler-Poisson systems. However, we can also view it as encoding the condition that our solution must expand; (48) implies that initially, u(t, η(t, x)) = x + ∂ t θ(log (1 + t), x), and thus if we have sufficient smallness on θ and its derivatives, the velocity of each particle stays close to the radial direction, which forces expansion. As discussed in the introduction, in the case of N ≥ 2, if the stars initially are far away from each other, and are then pushed further away, the gravitational interaction between two stars should be small. Thus the motion of each star should be close to that of constant velocity, as in [63], with the velocity essentially decomposing in to two main parts; one driving repulsion, and one driving expansion. For each κ = 1, 2, . . . , N , we define the error θ κ and repulsive velocity µ κ by the following relation This form of η κ encodes the fact that the particles should have a radial velocity to force expansion, and a repulsive velocity, to force the stars away from each other. As in [30,63], we define a new logarithmic timescale and under this transformation, we can write η κ as where ζ κ is defined for notational convenience. Under this ansatz, we will prove that solutions to (44)-(46) with small enough θ κ , as well as some more technical assumptions on µ κ , are global-in-time.
Remark 3.1. Note that the logarithmic timescale in (50) is used to study certain types of self similar solutions to evolution equations (see [17]), but in our case is for convenience.
We have some definitions to record. Throughout the paper, the identity matrix R 3 → R 3 will be denoted by I. Let Similarly to (39), we have for A [κ] and J κ : Using (50)-(51) to rewrite (44), we obtain and upon multiplying everything on the left hand side by e (1+ 3 α )τ , the system (44)- (46) becomes: Remark 3.2. The imbalance in exponential powers between the velocity term and the pressure term in (58) leads to a positive exponential in front of the velocity term in (59). This is exactly the nonlinear scaling mechanism, exhibited in [63], that will stabilise our solution and help us to prove that it is global-in-time.
For the gravitation potential term, note that in Eulerian coordinates, we have the Poisson equation (22) given by ∆φ = 4πρ, where ρ is the cumulative density κ ρ κ . Due to the compact support of each ρ κ , we can use the convolution formula for Poisson's equation and write φ explicitly as If we then apply the flow map η κ , for a fixed κ, to both sides of (62), we have, for i = 1, 2, 3, We can further rewrite G κ and where the last step uses integration by parts, noting thatw κ is 0 on ∂Ω. For I [κ,κ ′ ] , we have where the last line is because A [κ] = [∇ζ κ ] −1 . The G κ are the self-interaction terms and represent the part of the potential that encodes how the star is affected by its own gravity. The I [κ,κ ′ ] are the tidal terms which encode how different stars affect each other via their gravitational interaction. Note that in the case of N = 1 only the self-interaction terms are present, and these terms have already been studied in our context in [31]. The tidal terms are unique to the N ≥ 2 case, and understanding how to control these terms indirectly gives us information on how to configure our system initially, see Remark 6.4.

Initial Density Profiles
In this section we will fix our initial density profiles for each κ. First we define a collection of admissible profiles. Definition 3.3. For each κ ∈ {1, 2, . . . , N }, let W κ be the set of functions F : Ω κ → R with the following properties: • There exists a positive constant C > 0 such that for any x ∈ Ω κ where x → d(x, ∂Ω κ ) is the distance function to ∂Ω κ .
• The function given by is smooth on Ω κ . Now we fix an initial density profile for each κ = 1, 2, . . . , N by choosing a function from each W κ , which we will call W κ , and a constant δ > 0 such that Note that δ is significant as we can adjust the size of our initial densities by adjusting the value of δ. Throughout, δ ≪ 1, and we will specify explicit bounds on δ wherever necessary. This ansatz has already been used in [31,63] to find global-in-time solutions with small density. Under this definition ofρ κ , the speed of sound given byc 2 κ = γδW κ will satisfy the physical vacuum condition (17).

General Notation
For a function F : O → R, some domain O, the support of F is denoted supp F . For a real number λ, the ceiling function, denoted ⌈λ⌉, is the smallest integer M such that λ ≤ M . For two real numbers A and B, we say A B if there exists a positive constant C such that and for two real valued functions f and g, we say f g if f (x) g(x) holds pointwise. For two real-valued nonnegative functions f, g : O → R ≥0 , some domain O, we say f ∼ g if there exist positive constants c 1 and c 2 such that for all x ∈ O.
We also record the definition of the radial function on Ω = B 1 the unit ball: It is convenient to define shorthand for the distance function on Ω. Define

Derivatives
As we have seen above, rectangular derivatives will be denoted as ∂ i , for i in 1, 2, 3. In addition, we define various rectangular and ζ κ Lie derivatives that will be used throughout. The gradient, divergence, and curl on vector fields are given by for i, j = 1, 2, 3.
The ζ κ versions are given by In addition, we also need the matrix ζ κ curl, given by Recall (29); our reference domain is the closed unit ball in R 3 , Ω = B 1 . There exists a natural choice of spherical coordinates (r, ω, φ). An advantage of this choice of domain is that we can privilege the outward normal derivative, the direction in which the degeneracy of the problem occurs, due to the vacuum boundary condition (17). Accordingly we will, in essence, use ∂ r as the normal derivative, and ∂ ω , ∂ φ as the tangential derivatives. However we modify these derivatives by using linear combinations. These modifications allow for better commutation relations with the rectangular derivatives. Let the angular derivatives ✁ ∂ ij and radial derivative Λ be given by where the x i and ∂ j are rectangular, and i, j run through 1, 2, 3. On regions separated from the origin, we will use the following decomposition frequently: Remark 4.1. The coefficients of the derivatives we have defined go to 0 at the origin, which means we can only use them for estimates on a region separated from the origin. This can be dealt with using a partition of unity argument. Near the boundary we use these modified spherical derivatives, and on the interior, we are free to use rectangular derivatives as the degeneracy at the vacuum boundary is not an issue in this case.
Now, for m ∈ Z ≥0 , and n = (n 1 , n 2 , Although there are six non-zero ✁ ∂ derivatives to consider, ✁ ∂ ij = − ✁ ∂ ji , so (91) covers all cases. For such an n ∈ Z 3 ≥0 , |n| = n 1 + n 2 + n 3 . Similarly for rectangular derivatives we define, for k = (k 1 , We have the commutation relations between the modified spherical and rectangular derivatives, for i, j, k, m ∈ {1, 2, 3}, given by We also define commutators between the higher order differential operator defined in (91), and ∇: We can do the same thing with ∇ ζκ : There is no corresponding definition to (99) for ∇, as ∇ and ∇ k commute for all k ∈ Z 3 ≥0 . Note that (98) and (99) also define analogous objects for Curl ζκ and div ζκ as the former is ∇ ζκ − ∇ ζκ ⊺ , and the latter is Tr ∇ ζκ .

Energy Function
Following the strategy set out in Remark 4.1, define a smooth radial cutoff function χ on the closed unit ball such that In addition defineχ byχ Recalling the definition of α given in (27), we now define our energy spaces.
The norm of X b κ is given by with semi norm given by Remark 5.2. The powers ofW κ in the integrals involvingχ in (102) − (103) are not consistent with what we see in the definition of the energy spaces. However, on suppχ,W κ ∼ 1 so this discrepancy does not create issues.
We now define our higher order energy function and Curl ζκ energy function.
The Curl ζκ energy function of order b for (∂ τ θ κ , θ κ ) is given by Then the cumulative energy function and Curl ζκ energy function of order b are 6 Main Result and A Priori Assumptions

Main Result
In this section we will state the central theorem of this paper, global-in-time existence of the system (72)- (74). First we state the local-in-time theory for this system. Details on its construction are given in Appendix A.
Recall the definition of δ, the smallness parameter for the initial density profiles defined in (69), and the repulsive velocities of each star µ κ , defined in (49).
holds, and the bound holds. Then there exists a T > 0 such that we can find a unique solution Moreover, the function τ → S M (τ ) is continuous.
We introduce the Strong Separation Condition (SSC), that will be key in formulating our main result. holds.
Remark 6.3. We see that the Strong Separation Condition is a lower bound on any convex combination of the relative positions and velocities.
Remark 6.4. Note that the SSC implies the separation condition (108) automatically for large enough L. Moreover, we remark that (112) arises naturally from finding sufficient bounds on the tidal terms I [κ,κ ′ ] , defined in (63), to prove global-in-time existence of our solutions. Hence studying the tidal terms gives us crucial information on the initial geometry of our star configurations.

A Priori Assumptions
In this section we list all the assumptions we will make to prove Theorem 6.5. First we make explicit the assumption, implicit in (113), that ∇µ κ must be small in the X M κ norm for κ = 1, 2, . . . , N . Let 0 < ε 1 ≪ 1. We assume Let 0 < ε 2 ≪ 1. We specify the L for which the Strong Separation Condition (112) holds. Assume Now we state our a priori assumptions on the solution to (72)- (74). There exists a T > 0, depending on ε 2 , such that on the time interval [0, T ], the solution to (72)-(74) satisfies As a justification for why we can make our a priori assumptions, note that such a time interval [0, T ] where (117)-(121) hold must exist for small enough initial data by the local well-posedness theory set out in Theorem 6.1. Once we have used these assumptions to prove Theorem 6.5, we shall use the global-in-time energy bound (114) to improve our a priori assumptions, thereby justifying them via a continuity argument. Just like δ defined in (69), both ε 1 and ε 2 will be need to be small for our proofs. Throughout they will be ≪ 1, and where needed we will state any explicit bounds for them. Remark 6.6 (Asymptotic velocities). Note that assumption (113) means that each repuslive velocity µ κ is very close to a constant vectorμ κ . Thus we can write our initial velocities as As discussed in Section 3.1, x −x κ corresponds to expansion, andμ κ corresponds to repulsion. Moreover, the Lagrangian flow maps, and the Lagrangian velocities take the form Origin Initial displacement of centres from origin,x κ .
Expansion of stars for τ > 0.
Thus, any initial configuration of stars, with repulsive velocity equal to the initial displacement of their centre from the origin, launches a global-in-time solution as long as they are sufficiently separated so that (125) holds. Figure 1 represents a particularly symmetric example.

Estimates for the Gravitational Potential
In this section we obtain the estimates we need for the self interaction and tidal terms defined in (63), G κ and I [κ,κ ′ ] , to prove global-in-time existence of our solution. The estimates for G κ will be closely follow the methods used in [31], as they have to deal with the corresponding term in the One Body Euler-Poisson system. However, the tidal terms I [κ,κ ′ ] will clearly only appear in the case of two or more interacting bodies, and in the estimates for these terms we make crucial use of the Strong Separation Condition (117).

Tidal Term Estimates
Proof. We begin with the integrals on the left hand side of (126) that are localised on supp χ. For i = 1, 2, 3, we due to applying the Liebniz rule to Λ m ✁ ∂ n I i [κ,κ ′ ] , and the last term on the right hand side, I 3 , is every resulting term except for the cases when all the derivatives fall on either L (a, b, c, d) are constant coefficients resulting from the application of the Leibniz rule.
We will show the estimates for I 1 and I 2 , with I 3 following similarly. Since these tidal terms measure the gravitational interaction between the stars, it is natural that the separation of the bodies would influence our estimates, and indeed for both I 1 and I 2 , the Strong Separation Condition (117) is used. We begin with I 1 which is the simplest term, and then move on to I 2 which requires more care, especially when applying the Leibniz rule to Bound for I 1 .
We have that ζ κ (τ, where the first bound is due to the reverse triangle inequality, and the second is due to the Strong Separation Condition (117), with λ = e −τ , as well as the bound due to the a priori assumption (121). Thus For the second bound, we use (128) and the fact thatW κ is bounded on Ω for all κ to bound the z integral by a constant. Note that x + e −τx κ is smooth so we can bound derivatives of this term in L ∞ (Ω). We bound the θ κ term by the cumulative energy function S M . It remains to bound the µ κ on the right hand side.
When m + |n| = 0, we have that Here we have used (128), and the fact that χ andW κ are in L ∞ (Ω) for all κ ∈ {1, 2, . . . , N }. To bound ζ κ ′ L ∞ (Ω) we have Then we use (218) in Lemma D.4 for the θ κ ′ term, and (131) for the µ κ ′ term to obtain so to bound Λ m ✁ ∂ n (|ζ κ (x) − ζ κ ′ (z)| −3 ) effectively, it is enough to bound each of the terms in the sum separately, for every valid choice of p, a, b, c, d 12 , d 23 , and d 13 . Strictly speaking the correct expansion would include varying constants in front of every term in (136) depending on each index being summed over. They have all been set to 1 as they are not important when finding sufficient bounds for each term. We can write every separate term in (136) as where, once again, constants in front of each term in the sum have been set to 1. For any of the c, d 12 , d 23 or d 13 equal to 0, the product above is an empty product, equal to 1.
Finally we can expand the derivatives above to see that (137) can be written as a linear combination of terms of the form subject to the condition that Thus we see that Then note that we have For the first inequality, we first use (128) to bound |ζ κ −ζ κ ′ | −(3+2p) in L ∞ (Ω). We then apply L 2 (Ω) Cauchy-Schwartz to the z integral over Ω, and apply R 3 Cauchy-Schwartz to obtain To see that we can write the upper bound as on the second line in (142), notice that if g + |h| ≥ 1, then (143), from the resulting upper bound we can collect the p 1 terms that look like Λ g ✁ ∂ h ζ κ , and the 2p − p 1 terms that look like |ζ κ − ζ κ ′ |.

The second inequality in (142) is because Λ a1
has no z dependence and can be taken outside of the integral. The remaining integrand isW 2α κ ′ |ζ κ − ζ κ ′ | 2(2p−p1) |ζ κ ′ | 2 which can be bounded in L ∞ (Ω), using the assumptions onW κ ′ in Definition 3.3, as well as (131), (135), and the a priori assumption (118). Therefore we look to bound We can assume without loss of generality that 1 ≤ a 1 + |b 1 | ≤ · · · ≤ a p1 + |b p1 | ≤ m + |n|. We write For i = 1, . . . , p 1 − 1, we have a i + |b i | ≤ a p1 + |b p1 |, and therefore a i + |b i | ≤ (m + |n|)/2 ≤ M/2. Thus, we can apply . . , p 1 − 1, and bound these terms in L ∞ (Ω). We obtain The θ κ term on the right hand side can be bounded by the cumulative energy function S M . For the remaining term, we use (130) to bound the µ κ term, and note that we can boundW m− ai κ in L ∞ (Ω), as a i ≤ m. Thus we obtain where the last bound is due to the a priori assumption (118).
The bounds (133) and (147) cover all possible cases for I 2 , thus we have as required.
Bound for I 3 .
Similarly to I 2 , to bound any of the terms in the sum that forms I 3 it is sufficient to bound all terms of the form subject to the condition that a + a i = m, and |b| + |b i | = |n|. This requires exactly the same strategy as (146), and so we immediately obtain The bounds in (129), (148), and (150), along with analogous estimates on suppχ give the proposition.
Now we move on to estimates for the gravitational potentials acting on each separate body coming from their own mass.

Self Interaction Term Estimates
In this section we estimate the term coming from the gravitational effect of a body in this system on itself, G κ . Let us recall that we can write The key estimate for G κ is laid out in the following proposition.
We begin with estimates for the integrals near the boundary, on supp χ. The strategy follows the strategy in [31]. We first bound tangential derivatives of G , and use these, as well as the Poisson equation ∆φ = ρ, to close estimates in the normal direction. The esimtates for tangential derivatives of G κ near the boundary are stated in the following proposition: Proof. The proof relies on the following lemma, which was first used by Gu and Lei [25] in the context of the spatial domain being T 2 × R rather than the closed unit ball. It was then adapted in to our setting by Hadžić and Jang [31]. For a full proof of the proposition, see [31], whose methods apply directly to our case.
Lemma 7.4. Let ✁ ∂ x and ✁ ∂ z denote angular derivatives in the x and z variables respetively. Let h 1 : R 3 → R, and h 2 : for some constants c ν ′ .
Proof. The proof is a combination of the identity ✁ ∂ x = ✁ ∂ x + ✁ ∂ z − ✁ ∂ z , as well as integration by parts in the z variable.
We now use Proposition 7.3 to prove Proposition 7.2. Once again, the proof strategy very closely follows that of the one in [31].
Proof of Proposition 7.2. Recall that we are looking to prove We first concentrate on the interior bounds, localised on suppχ. Recall the definition of χ andχ as smooth radial cutoff functions in (100) and (101). We know that suppχ = B 3 4 , the ball of radius 3/4 around 0. Define another smooth radial functionχ ν such that such that 3 4 + ν ∈ (3/4, 1). Correspondingly we define To bound the integrals on suppχ we have For the first term on the right hand side above, we can use the methods we used in Proposition 7.1 as for x ∈ suppχ, z ∈ supp χ ν , using the mean value theorem and (119). For the second term we can use the same methods as Proposition 7.3, noting thatχ ν and its derivatives are 0 on ∂Ω, so there are no boundary terms when we integrate by parts to get a obtain identity to Lemma 7.4. Therefore, we have Now we concentrate on the supp χ bounds. As in [31], we use induction on the number of radial derivatives. For now we fix a κ. Our induction hypothesis is where C(a) is a constant depending only on the number of radial derivatives and our a priori assumptions (116)-(120), and ε 1 , ε 2 are the constants mentioned in the a priori assumptions, Section 6.2. The case where m = 0 has been dealt with in Proposition 7.3. Now assume it holds for some 0 < a < M . Following [31], we use the Poisson equation (22) to write ΛG κ in a form amenable to the necessary estimates.
Recall that (22) gives us the relation ∆φ = 4πρ, and we know from (15) that ρ is the cumulative density defined by ρ = κ ρ κ , so, where ρ, {ρ κ } are thought of as functions on R 3 with compact support, we can define φ κ to be the solution to on R 3 , with φ κ → 0 as |x| → ∞ as the boundary condition. We know that φ κ = Φ * ρ κ , where Φ is the fundamental solution, so in Lagrangian coordinates we have where we recall the definition of G κ from (63). We take the Lagrangian divergence of this relationship, and comparing to the pullback via η κ of (161), obtain Then, using (43), (47), (51), and (54), we have Moreover, since G κ can be written as the Lagrangian gradient of some potential function, we immediately have Now we define R κ : With this definition of R κ we have We also have Thus Then, using (170) and the commutativity of Λ and ✁ ∂, we obtain The first bound is due to (170) and the second is due to the induction hypothesis (160). It is left to estimate the div G κ and Curl G κ terms. For this we use (167) and (168). First we have If a 1 + |b 1 | = 0 then bothW α , and J −1 κ can bounded in L ∞ (Ω), due to the assumption that ∇W κ ∈ X M κ and the a priori assumption (120) respectively. Otherwise, we first use (208) (resp. (209)) from Lemma (B.2) to bound the terms involving L 2 (resp. L ∞ ) norms of derivatives of J −1 κ . Then for theW α κ terms, we can use ∇W κ ∈ X M κ (resp. (220) in Lemma D.4) to bound the L 2 (resp. L ∞ ) terms after using the Leibniz rule on Λ ai ✁ ∂ b i W α κ . Note here that our range of γ is crucial in keeping the norms of the derivatives ofW α κ finite. Next we have The terms involving L 2 (resp. L ∞ ) norms of derivatives of A [κ] can be bounded using (208) (resp. (209)) in Lemma B.2. For the G κ terms, we have for i = 1, 2: In (174) we use (90), and in (175) we use (90) and (219) in Lemma D.4. The R κ terms can be bounded analogously to (174) and (175). Thus Then, from the identities for div G κ and Curl G κ given in (167) and (168), as well as (171)-(176), we have for some constant C(a + 1). This proves the induction hypothesis.
Clearly the sum all of the C(a) can be bounded above by an absolute constant that we call C 1 , and so summing these bounds for all a, we obtain Combining (178) and (159), we have, for some constant Recall in Section 6.2 we defined ε 1 and ε 2 as small constants that could be shrunk if necessary. Upon doing this, and recalling the definition of X M κ norm from Definition 5.1, we have Summing over κ for (180) gives (155), and completes the proof of Proposition 7.2.

Potential Term Estimates
We end this section with the main estimate needed to bound all potential terms when we perform our energy estimates.
Proof. The proof is follows from the identity (63), and Propositions 7.1 and 7.2.
This is structurally very similar to the Curl ζκ equation obtained in [63], and indeed the subsequent estimates follow an analogous strategy. One important difference is that in our case, derivatives of A [κ] and J κ produce derivatives of µ κ as well as derivatives of θ κ . However, these terms can be dealt with using the assumption on ∇µ κ , (116). For example, we can bound terms like so: where we use (116) to bound the µ κ norm above by a constant. Hence, these extra terms are readily bounded without any further techincal difficulties. Therefore, using the methods in [63], we obtain the following theorem.

Energy Estimates
In this section, we prove the main energy identity, and subsequent energy estimates that will form the bulk of the proof of Theorem 6.5. First we need some definitions.
Remark 9.2. Just as in [30,63], this damping functional does not play an essential role in our analysis. We only require they have correct sign so that our energy estimates are sufficient. This is guaranteed by the fact that β = 3(γ − 1) ≤ 3/2, implied by our range of γ.
Similarly to [30], we define a truncated-in-time higher order energy function which we use to prove our main theorem.
The cumulative truncated energy function is then given by Note that we have the identity and, for τ 1 ≤ τ 2 , the inequalities With the damping functional from Definition 9.1, we can state our main energy identity.
with R κ (m, n) and R κ (k) being remainder terms that we can bound effectively.
Proof. The proof of Theorem 9.4 can be adapted from the methods used in [63]. For each fixed κ, on supp χ we first divide byW α κ in (72), and then act on the result withW α+m κ Λ m ✁ ∂ n , some 0 ≤ m + |n| ≤ M , resulting in the identity From here, we take the L 2 (Ω) inner product of this equality with χΛ m ✁ ∂ n ∂ τ θ κ . We can follow a similar procedure on suppχ withW α κ ∇ k andχ∇ k ∂ τ θ κ . From here we can extract the left hand side in (191), and move everything else to the right hand side. Doing this for each κ = 1, 2, . . . , N and then summing gives us the energy identity (191), once we have labelled the remainders of this procedure from the pressure term R κ (m, n) and R κ (k) respectively.
for some constants C 1 , . . . , C 5 ∈ [1, ∞), and G : 1 δ e βs ∂ τ θ κ (s) For each κ most of the terms coming from R κ (m, n) and R κ (k) will be analogous to the ones obtained in [63]. However, just as with the curl estimates in Section 8, we will also have remainder terms that look, for example, like a weighted L 2 (Ω) inner product of Λ m ✁ ∂ n µ κ and Λ m ✁ ∂ n ∂ τ θ κ , coming from the derivatives of A [κ] and J −1/α κ . These terms can be bounded, for example like where we use (116) to bound ∇µ κ X M κ by a constant.
Let T be such that for some constantC, whose existence is guaranteed by the Local Well-Posedness theory set out in Theorem 6.1. Let C * be defined by with C i , i = 1, 2, 3, 4, 5 given in Theorem 9.5. Note thatC < C * , so given we have T ≤ T * . Now letting τ 1 = T /2, Theorem 9.5 tells us that for any τ ∈ [ T 2 , T * ), we have Let δ satisfy Then we have Combining (203) with (197), and using (190) we obtain for all τ ∈ [0, T * ). Shrinking δ further if necessary, we also improve our a priori assumptions (118)-(121). For example, for all τ ∈ [0, T * ), for small enough δ. Similarly for J κ , and θ κ . Then, by continuity of S M (τ ) as a function of τ , we must have that T * = ∞. Therefore, the bound (114) follows. It is left to prove (115). Fix κ ∈ {1, 2, . . . , N }. Let τ 1 > τ 2 . For any m + |n| ≤ M , we have the estimate The bound follows from Cauchy-Schwarz in τ and the global energy estimate (114). A similar bound holds for any |k| ≤ N on suppχ. Thus θ κ (τ n ) is Cauchy, for any strictly increasing sequence τ n . As X M κ is a Banach space, lim τ →∞ θ κ (τ ) exists in X M κ ; we call this limit θ Proof. The full proof of (208) can be adapted very easily from [63]. The main technical difficulties are for the integrals localised near the boundary, on supp χ, with regards to ensuring each term resulting from the application of the Leibniz rule to either Λ m ✁ ∂ n (A [κ]) or Λ m ✁ ∂ n J −1 κ appears with a high enough power ofW κ to be bounded in either the X M κ or Y M κ (∇ ζκ ) norms. For the L ∞ estimates in (209), we first concentrate on the A [κ] term. Note that (56) gives and similarly for ΛA [κ]. We apply (210) repeatedly to obtain where l a+b , l νj +λ i are nonnegative integers. By (119), we have A [κ] L ∞ (Ω) 1. It is left to bound terms of the , for any ν ∈ Z ≥0 , λ ∈ Z 3 ≥0 such that 1 ≤ ν + |λ| ≤ m + |n|. First we have For the θ κ term, we write Λ ν ✁ ∂ λ ∇ as a linear combination of operators of the form ∇Λ ν ′ ✁ ∂ λ ′ , ν ′ + |λ ′ | ≤ ν + |λ| with smooth coefficients using (88), (89), and (90). Finally, employing (220) and (219) from Lemma D.4 to the θ κ and µ κ terms respectively gives us The proof of the same bound for the J κ term is similar, where we utilise the differentiation formula (57) instead.
with norm given by The definition of H s (O) can be extended to s ∈ R ≥0 by interpolation.
Definition D.2. For a bounded domain O ⊂ R 3 , α > 0 and b ∈ Z ≥0 , define the weighted Sobolev space H α,b by with norm given by where d = d(x, ∂O) is the distance function to the boundary on O.
Given these definitions, we have the following embedding. The proof of Lemma D.3 is standard and can be found in [43]. Finally, we state the Hardy-Sobolev bounds for L ∞ -norm in terms of our energy norms. Recall the spaces X b κ and Y b κ (∇ ζκ ) with their associated norm and semi-norm respectively, given in Definition 5.1. We also recall the definitions of χ andχ given in (100) and (101).