Central configurations in planar $N$-body problem for $N=5,6,7$ with equal masses

We give a computer assisted proof of the full listing of central configuration for $n$-body problem for Newtonian potential on the plane for $n=5,6,7$ with equal masses. We show all these central configurations have a reflective symmetry with respect to some line. For $n=8,9,10$ we establish the existence of central configurations without any reflectional symmetry. The proof is based on the interval arithmetics.


Introduction
In this paper, for n = 4, 5, 6, 7 we consider two questions: -finding all central configurations (CCs) in Newtonian n-body problem on the plane with equal masses and -showing that each central configuration has a line of reflection symmetry.

The state of the art
The listing (apparently full for n 9 ) of central configurations with equal masses was given by Ferrario in unpublished notes [F02] for n ∈ {3, . . . , 10} and by Moeckel [Moe89] for n 8. For n = 4 it was shown by Albouy that all CCs have some reflectional symmetry [A95] and later in [A96] with computer assistance the full list of central configurations was given. From numerical simulations (see for example [Moe89,F02]) it is apparent that all CCs with equal masses posses some symmetry for n = 5, 6, 7. Moeckel [Moe89] have found numerically some CCs without any symmetry for n = 8. Also for n = 9 Simó [Si] has found 2 families, non-equivalent, and without any symmetry. Some CCs without symmetry for n = 10 can be found also in [F02]. The investigations of central configurations for equal masses is a subcase of more general problem of central configurations with arbitrary positive masses. The general conjecture of finiteness of central configurations (relative equilibria) in the n-body problem is stated in [W41] and appears as the sixth problem of Smale's eighteen problems for the 21st century [Sm98]. There are many works on the existence of some particular central configurations. Here we discuss only those papers which aim to more general statement about all CCs. The two most important works are [HM05] and [AK12]. In [HM05] the finiteness of CC for n = 4 for any system of positive masses was proved with computer assistance. In [AK12] for n = 4 problem the finiteness of CC was proven without computer assistance. In the same paper, the finiteness for n = 5 was proven for arbitrary positive masses, except perhaps if the 5-tuple of positive masses belongs to a given codimension 2 subvariety of the mass space. It is interesting to notice that the equal masses case treated in our paper does not belong to this subvariety. For the spatial 5-body problem Moeckel in [Moe01] established the generic finiteness of Dziobek's CCs (CCs which are not planar). A computer assisted work by Hampton and Jensen [HJ11] strengthens this result by giving an explicit list of conditions for exceptional values of masses. A common feature of these works is that they give a quite poor estimate for the maximum number of central configurations. In this context it is worth to mention the work of Simo, based on extensive numerical studies. In [Si78] he gives the number of CCs for all possible masses for n = 4.
In [LS09] the spatial 5-body problem with equal masses was considered. A complete classification of the isolated central configurations of the 5-body problem was given (this includes also a planar isolated CCs). The approach has a numerical component, hence it cannot be claimed to be fully rigorous. Also the proof does not exclude the possibility that a higher dimensional set of solutions exists. On the other hand the existence of identified isolated CC, has been proven using the Krawczyk's operator, i.e. a tool from interval arithmetic we also use. Kotsireas (see [Ks00] and references given there) considers the 5-body problem with equal masses. He gives computer assisted proof of a full list of all such configurations and shows that each of them posses some reflectional symmetry.
The above mentioned works study the polynomial equations derived from the equations for CC using the (real or complex) algebraic geometry tools. In contrast, we take a different approach: we use standard interval arithmetic tools, hence in principle we can treat also other potentials which cannot be reduced to polynomial equations.

The main results
Theorem 1 There exist only a finite number of various types of CCs, for n = 5, 6, 7 the planar n-body Newtonian problem with equal masses. They are listed in Section A.1. Any CC can be obtained from one of them by suitable composition of translation, rescaling, rotation, reflection and permutation of bodies. Moreover, each of these central configurations has some reflectional symmetry.
Theorem 2 For n = 8, 9, 10 in the planar n-body Newtonian problem with equal masses, there exist CCs without any line of reflectional symmetry. They are listed in Section A.2.
In the case of equal masses one can consider equivalences in two different ways: either one passes from a solution to another one by rotation (scaling is already taken into account) or one can also add permutations and reflections. For instance, for 4 bodies in the first criterion of equivalence there are 50 classes (see numerical work by [Si78]), while in the second only 4 classes. In this paper we use this second criterion for equivalence.
Let us give some comments about our method. This is basically a brute force approach using standard interval arithmetic tools. Looking for CCs we consider the whole configuration space, and it is surprising that the most demanding part is to exclude the possibility of the existence of CC in a given box. Once we are 'very close' to an isolated CC, it is relatively easy to establish its existence and local uniqueness using the Krawczyk's operator [K69]. The additional difficulty is that the potential contains singularity, which introduces some non-compactness in the domain to be covered. Consequently our algorithm, which is more or less a binary search algorithm, scales poorly with n -this is the so called dimensionality curse (see [TWW88]), which means that the complexity of our algorithm grows exponentially with n. For example, assume that we can decide if a box contains a body belonging to some CC, only when it is less then 10 −2 in each direction. Then adding a new body in [−1, 1] × [−1, 1] multiplies the number of boxes to be examined by (2/10 −2 ) 2 = 4 · 10 4 . Hence we do not include a rigorous listing of CCs for n = 8. Note that for n = 5 the computations were done in 24 seconds, for n = 6 it took about one hour to get the result, while for n = 7 we needed almost a hundred hours.
For any CC from the listing in [Moe89] or [F02] for n 8 we have found no difficulty proving its existence and local uniqueness. In particular we confirmed the existence of non-symmetric planar CCs for n = 8, 9, 10 (see Theorem 2).
The paper is organized as follows. In Section 2 we recall the equations for the central configurations and their basic properties. In Section 3 we derive several a priori bounds for CC, so that we obtain a compact domain for the search algorithm we implemented. In Section 4 we discuss various tests which are used to show that a given box does not contain any CC. In Section 5 we derive a reduced set of equations for CC. This is necessary because to apply Krawczyk's method we need to ensure that the system of equations does not contain any degeneracies, which are due to symmetries of the original system of equations for CCs. In Section 6 we give assumptions and basic ideas concerning the computer assisted proofs of main Theorems 1 and 2. We explain the Krawczyk's method. Details of the algorithm are described in Section7. In Section 8 we present an attempt to minimize dependency problem in interval arithmetic in evaluation of gravitational force. In the Appendix we give an output of the program establishing Theorems 1 and 2 (also for n = 3, 4) and pictures of all CCs found.

Equations for central configurations
Let q i ∈ R d , i = 1, . . . , n and d 1 (the physically interesting cases are d = 1, 2, 3), where q i is a position of i-th body with mass m i ∈ R + . Let us set Central configurations are solutions of the following system of equations (see [Moe]): where λ ∈ R is a constant, c = ( n i=1 m i q i ) /M is center of mass, r ij = r ji = |q i − q j | is the Euclidean distance between i-th and j-th bodies and (−f i ) is the force which acts on i-th body resulting from the gravitational pull of other bodies. The system of equations (2) has the same symmetries as the n-body problem. It is invariant with respect to group of isometries of R d and the scaling of variables.
In the planar case if we consider the bodies in a rotating system (with the center of mass at the origin) with constant angular velocity ω = √ λ, the physical meaning of (2) is obvious: the gravitational attraction is compensated by the centrifugal force and the central configurations are fixed points in the rotating frame. For more on the origin and physical meaning of equations (2) the reader should consult [Moe] and the references given there.
The system (2) has dn equations and dn + 1 unknowns: q i ∈ R d for i = 1, . . . , n and λ ∈ R + . The system has a O(d) and scaling symmetry (with respect to q i 's and m i 's). If we demand that c = 0 (which is obtained by a suitable translation) and λ = 1 (which can be obtained by rescaling q i 's or m i 's) we obtain the equations (compare [Moe, Moe14, AK12]) 1 m i f i (q 1 , . . . , q n ), i = 1, . . . , n.
It is easy to see that if (3) is satisfied, then c = 0 (see Sec. 2.1) and (2) also holds for λ = 1. A q = (q 1 , . . . , q n ) ∈ R d n is called a configuration. If q satisfies (3) then it is called a normalized central configuration (abbreviated as CC). For the future use we introduce the function F : Then the system (3) becomes F (q 1 , . . . , q n ) = 0.

Some identities and conservation laws
It is well know that for any (q 1 , q 2 , q 3 , . . . , where v ∧ w is the exterior product of vectors, the result being an element of exterior algebra. If d = 2, 3 it can be interpreted as the vector product of v and w in dimension 3. The identities (6) and (7) are easy consequences of the third Newton's law (the action equals reaction) and the requirement that the mutual forces between bodies are in direction of the other body. But (6) and (7) can be seen also as the consequences of the symmetries of Newtonian nbody problem. According to Noether's Theorem, by the translational symmetry we have a conservation of momentum, which is equivalent to (6), while the rotational symmetry implies the conservation of angular momentum, which is implied by (7).
Note that the components of v ∧ w are given by determinants. In any dimension in the presence of the rotational symmetry, for any direction of rotation identified by v 1 ∧ v 2 ( v 1 and v 2 are perpendicular unit vectors) the following quantity must be zero (as a consequence of the Noether Theorem and the invariance with respect to the rotation in the plane Consider system (3). After multiplication of i-th equation by m i and addition of all equations using (6) we obtain (or rather recover) the center of mass equation We can take the equations for n-th body and replace it with (9) to obtain an equivalent system.
Later in Section 5 we will use (7) to defined a reduced system of equations for CCs which will not have the degeneracies present in (3).

Moment of inertia of central configurations
The important role of the moment of inertia in the investigation of central configurations is well known. In our context it plays a crucial role in stating some a priori bounds for central configurations. We present, with proofs, some well known results on moment of inertia taken from the notes by Moeckel [Moe14] and the paper of Albouy and Kaloshin [AK12].
Definition 1 For a configuration q let the moment of inertia I(q) and the potential function U (q) be given by Lemma 3 Assume that i m i q i = 0 and M = 1. Then Proof: Let us denote q i,j = q i − q j . Since we have Lemma 4 If q ∈ R d n is a (normalized) central configuration, then Proof: We take the scalar product of i-th equation in (3) by m i q i and add these equations to obtain

A priori bounds for central configurations
From the point of view of CAP (computer assisted proofs) in the problem of finding and counting all CCs the issue of compactness of the search domain is fundamental. The lack of compactness arises in two cases: -two or more bodies might be arbitrary close to a collision, -some bodies might be arbitrary far from the origin.
The goal of this section is to deal with these issues. We will give a priori bounds (depending on m i 's) on the minimal distance of the closest bodies and for the size of the central configuration.

Lower bound on the distances
It is well known that central configurations are away from the collision set (see [Sh70] or [Moe14,Prop. 15]). However in these works no quantitative statement directly applicable to system (3) has been given. Here we develop explicit a priori bounds. The main idea is to use I(q) = U (q) (see Lemma 4) to show that some term(s) m i m j /r ij entering U (q) dominate and cannot be balanced, when bodies are very close. Observe that using I(q) = U (q) and positivity of all terms entering U (q) allows us to escape the discussion of large terms on the rhs in the system (3), which might cancel out or not etc. This is not the case in the framework of Albouy and Kaloshin [AK12], where complex configurations and even complex masses have been considered, hence the positivity of I and U is lost.
Here we present the simplest result; later we'll give more accurate criteria.
Lemma 5 Assume that a (normalized) central configuration q ∈ R d n satisfies |q i | R for i = 1, . . . , n. Then I(q) M R 2 .
Proof: Since for any 1 i n : Theorem 6 Assume that a (normalized) central configuration q ∈ R d n satisfies |q i | R for i = 1, . . . , n. Then Proof: From Lemma 4 and Lemma 5, for any 1 i < j n, we have hence (15).
Below we establish a lower bound on the radius of ball centered at 0 and containing a central configuration in the case of equal masses.
Theorem 7 Assume that all masses are equal m and q is a normalized central configuration such that |q i | R. Then Proof: Since for any 1 i n : |q i | R, thus r ij 2R and we obtain the following bound Hence from Lemma 4 and Lemma 5 we obtain If M = 1, then lim n→∞ 3 n−1 4n = 4 −1/3 ≈ 0.629961. This estimate appears to be reasonably good, as shown in Table 1.  Table 1: The size of a minimal ball containing all normalized central configurations with M = 1 for several n's for equal masses case. The minimum R is realized for regular n-gon.
In the next theorem we do not assume that all masses are equal.
Theorem 8 Assume M = 1 and q is a normalized central configuration. Then there exists a pair i = j such that r ij 1.
Proof: For the proof by contradiction assume that r ij < 1 for all pairs of bodies. Hence we have r 2 ij < 1/r ij for all pairs. From Lemma 3 it follows From Lemma 4 it follows that q is not a central configuration.
Observe that the above estimate is optimal, because it is realized for the equilateral triangle for n = 3 and a tetrahedron (non planar CC) for n = 4. From the above theorem we obtain the following lower bound for the size of a central configuration. Contrary to Theorem 7 we do not assume that all masses are equal.
Theorem 9 Assume that M = 1 and q is a normalized central configuration and |q i | R for i = 1, . . . , n. Then R 1/2.
Proof: For the proof by contradiction assume that R < 1/2. Then for all pairs r ij < 1. From Theorem 8 it follows that q is not central configuration.

The upper bound on the size of central configuration
The goal of this section is to give the upper bounds for the size of the central configuration. This time we exploit the fact that if the forces are bounded, then large q i 's on the left hand side of the system (3) cannot be balanced. The obvious difficulty with the realization of this idea is: we can have a group of bodies with large norms which are close to each other in the central configuration, which produce large terms on rhs of the system (3). To overcome this we consider clusters of points far from the origin and the resulting force on it. In such situation mutual interactions between bodies in the cluster cancel out.
Lemma 10 Assume q ∈ R d n is a central normalized configuration. Let R = |q i 0 | = max i=1,...,n |q i |. Then for all ε ∈ (0, R/(n − 1)) holds Proof: For simplicity let's assume that d = 2 and q i = (x i , y i ). Let us fix any ε ∈ (0, R/(n − 1)). After a suitable rotation of coordinate system we can assume that Let C be a minimal subset (cluster) of indices of bodies satisfying the following conditions The cluster C can be constructed as follows: We start with i 0 ∈ C. Then we add all bodies which are not farther than ε from the bodies already in C. We repeat this until the set C stabilizes, which will happen after at most n − 1 steps. From assumption about ε and R it follows that Observe that (20) implies that C = {1, . . . , n}. Indeed (20) and (19) imply that x i > 0 for all i ∈ C. This and the center of mass condition (9) implies that C cannot contain all bodies. This implies that the process of building C must stop after at most n − 2 steps. Therefore we obtained a cluster C with the following properties q i0 R r Figure 1: R > (n − 1)ε and r = (n − 2)ε; the darkened area is the region where all the bodies from cluster are located.
We have For n 10 and the equal masses the above estimate appears to be an overestimation. The largest possible size of CC was slightly above 1 and is considerably smaller that one established in the above theorem.

Exclusion tests for CCs
Assume that we have an interval set D (i.e. a box, a product of intervals) on which we would like to exclude a solution. We do not assume that D ⊂ dom F (see (5)) and this is an important point. The a priori estimates discussed in Section 3 allow to exclude D iff there is no point in D which satisfies the obtained bounds.
In the following we discuss other exclusion tests.

Checking for zeros
One obvious test is to check whether in the interval evaluation of F (D) (see (5)) at least one of the component does not contain zero. Observe that a partial tests of this type are possible even if D admits the collision, i.e. formally it is not contained in dom F , but we can do this verification whenever D ⊂ dom F i . The test takes the following form: 0. given box D in the configuration space, such that D ⊂ dom F i (i.e. the i-th does not have a collision with other bodies). Let D i = {q i | q ∈ D}.
1. compute f i (D) the interval enclosure of f i (q) for q ∈ D, (3)), then D does not contain any normalized central configuration Observe that point 2 allow us exclude the box D, but if this is impossible then we can replace D by D ref obtained in point 3, as any CC contained in D must belong to D ref . This is of several places in our algorithm, where we attempt to do better than doing naive binary subdivision in order to relieve the curse of dimensionality.

The cluster of bodies -checking for zeros test
In the case of colliding bodies or a cluster of close bodies, say with indices i = 1, . . . , s, after adding first s equations multiplied by m i we obtain Therefore for the cluster of bodies C ⊂ {1, . . . , n} we will check whether where the expression on the rhs of (29) is evaluated in the interval arithmetic on D.
If it is not satisfied then we conclude that D does not contain any CC. Again let us stress that the set D might contain some collisions, and this test is still applicable.

The cluster of bodies -test based on moment of inertia and potential
We will exploit I(q) = U (q) for CC (compare Lemma 4), but our focus will be on the subsets (clusters) of bodies. Let us fix C ⊂ {1, 2, . . . , n} and Z ⊂ R d n . Let us define In the case when C = {1, . . . , n} we set F C,Z = 0. The important point is that we can compute the infimum in U C,Z even if the set Z contains collisions. It makes sense to take as C a cluster of close points (containing possible collisions and near collisions), so that there is no collision between bodies in C and its complement. In such case F C,Z will be finite.
We have the following criterion for nonexistence of CC in Z: Lemma 12 Assume that U C,Z , I C,Z , F C,Z ∈ R and then there is NO central configuration in Z.
Proof: Without any loss of the generality we can assume that C = {1, . . . , s} for some 1 s n. Consider system (3). We multiply i-th equation by m i q i and we add first s equations. We obtain (compare the proof of Lemma 4) Let us stress that (30) must hold for any central configuration. Now for q ∈ Z holds Hence (30) is not satisfied. Therefore we do not have any central configuration in Z.
Observe that if C = {1, . . . , n} the above lemma is reduced to checking whether inf q∈Z U (q) > sup q∈Z I(q).

Where it should be helpful?
Assume that we have a cluster of bodies C with small moment of inertia and big potential energy and that there are bodies away from zero, thus it is possible to have I(q) = U (q). Then the above criterion should be effective to rule out such CCs. Definition 2 We will say that a normalized central configuration q = (q 1 , . . . , q n ) is nondegenerate if the rank of DF (q) is equal to dn − dim SO(d). Otherwise the configuration will be called degenerate.
The idea of the above notion of degeneracy is to allow only for the degeneracy related to the rotational symmetry of the problem, because by setting λ = 1 in (2) and keeping the masses fixed we removed the scaling symmetry.

The reduced system on the plane
We consider a planar case here, i.e. d = 2. The fact that the system of equations (3) is degenerate (each solution give rise to a circle of solutions) make this system not amenable for the use of standard interval arithmetic methods (see for example the Krawczyk operator discussed in Section 6.3)) to rigorously count all possible solutions. We need to kill the SO(2)symmetry and then hope that all solutions will be non-degenerate. In this section we show how to reduce the system (3) to an equivalent system amenable to the interval analysis tools.
Let us fix k ∈ {1, . . . , n − 1} and consider the following set of equations where f i = (f i,x , f i,y ). Observe that the system (32-34) coincides with the system (10-11) with the equation for y k dropped. The next theorem addresses the question: whether from the reduced system (32-34) we obtain the solution of (3)?
Theorem 13 Assume d = 2. If (q 1 , . . . , q n ) is a solution of the reduced system (32-34) and x k = x n , then it is a normalized central configuration, i.e. it satisfies (3).
Proof: Let q i be as in our assumptions. We need to show that m k y k = f k,y . From (34), (6) and (7) it follows that Since from (32) Let us take the look at n−1 i,j=1,i =j m j f i ∧ q j . We have from (32) hence again from (32) it follows that From the above and (35) we obtain From (33) we have Now if x k − x n = 0, then m k y k = f k,y .
The system (32-33) contains 2(n − 1) − 1 equations in 2(n − 1) variables and has O(2) symmetry (i.e. rotations around origin and reflection symmetries with respect to lines passing through the origin map solutions of this system into itself). In order to obtain a system with the same number of equations and variables we can impose additional condition leading to the removal of y k variable, so that the rotational symmetry will be broken. Obviously in the above setting we could drop the equation for x k and we will obtain an analogous result.
We can think of a general reduced system as follows: • we fix some hyperplane H, in the reduced (by the center of mass condition) configuration space R 2(n−1) , so that H is transversal to the action SO(2) and k is such that v k ∈ {x k , y k } can be computed in terms of other variables. This will induce an embedding, • in the system (10-11) we remove the equation for v k . Then the reduced system can be written as where R k is a projection which removes v k variable in the image.
The system (32-34) supplemented by substitution y k = y k (. . . ) is an example of (36). We present the following obvious result Theorem 14 Assume for simplicity that H is given by Assume that q is a solution of reduced system (36) with a substitution (37), such that x k = x n . Then q is a normalized central configuration. If q is non-degenerate solution of the reduced system, then this is non-degenerate solution of (3).
Proof: The first part is obvious in view of Theorem 13 and condition x k = x n implies that it is a central configuration. The maximum rank in the reduced system gives the non-degeneracy of the configuration in the sense of Definition 2.

Following [AK12] we have tried two possibilities
• we set k = 2 and we eliminate variable y 2 by setting • we set k = n − 1 and we eliminate variable y n−1 by setting Observe that for both normalizations defined above for any CC q there is a rotation R such that Rq satisfies this normalization. Hence we can safely impose any of those conditions without losing any CC. In both cases we will need First consider (38). Condition (40) does not hold for some CC in the case of equal masses. For example, for n = 4 and an equilateral rectangle, such that x 2 = x 4 we obtained numerically (and also symbolically using Mathematica) that the jacobian matrix for the reduced system has a zero eigenvalue. Hence the solution is degenerate for the reduced system. Now, consider the condition (39). If we setup our computations so that q n−1 body maximizes the distance from the origin for all bodies, then we have (40) satisfied, otherwise q n will be further from zero. This observation does not prove that if q is a non-degenerate CC in the sense of Definition 2, then it is also a non-degenerate solution of the reduced system, but this appears to happen in our rigorous computation of central configurations so far.

On the computer assisted proof
We normalized masses so that M = i m i = 1. In this section we index bodies from 0 to n − 1 to be in the agreement with our program. In the sequel we study the following reduced system where we set 6.1 Equal mass case, the reduction of the configuration space for CCs In the case of equal mass, after a suitable rotation and permutation of the bodies, we can assume that Condition (45) guarantees that x n−2 = x n−1 , hence the solution of a reduced system (41-42) is CC. In view of symmetry and Lemma 9 we impose some more conditions After suitable permutation of bodies and a reflection with respect to 0X-axis it is easy to see that each CC has its equivalent in the set of the configurations satisfying the following conditions • q n−2 = (R, 0) is the furthermost body from the origin • q 0 is the leftmost with non-negative y-coordinate • q 1 has the smallest y coordinate • all other bodies have their x coordinates in the order of increasing/decreasing indices.
This, combined with Lemma 10, shows that it is enough to consider the following set in which we look for the central configuration x 2 x 3 · · · x n−3 x n−1 .
We call this order increasing due to the requirement (53). In the computation we use analogous decreasing ordering in which we state the opposite, i.e.
For now, we do not know why it is better to use the decreasing order, but the difference is significant (see Table 2).
n increasing decreasing 4 0.027170 0.026153 5 3.369141 2.615399 6 1103.138988 924.083085 Table 2: Times of asynchronous computations in minutes for different orderings (the computations were carried out on the computer Intel Core i7-5500U CPU @ 2.40GHz 4 with 8GB RAM; a single thread was used).

Outline of the approach
In the algorithm we look for all zeros of the reduced system (41,42), which under assumption x n−1 = x n−2 is equivalent to (3). For our algorithm proving an existence of locally unique solution in some cube is as important as proving that in a given box there is no solution.
For proving of the existence of the locally unique solution we use the Krawczyk method applied to the system (41,42). To rule out the existence of solution we use the exclusion tests discussed in Section 4 and also the Krawczyk operator.
The symmetry of CCs is established by proving the uniqueness in a suitable symmetric box (see Section 7.2).

The Krawczyk operator
The Krawczyk operator[Al94, K69, N90] is an interval analysis tool to establish the existence of unique zero for the system of n nonlinear equations in n variables. Below we briefly explain how the Krawczyk operator is derived, as it appears mysterious and little known outside the interval arithmetic community.

Motivation, heuristic derivation
Let F : R n → R n be a C 1 -function. We would like to solve the equation We begin by explaining the basic idea of the Krawczyk method. The Newton method is given by It is well known that if F (x * ) = 0 and dF (x * ) is nonsingular, then x * is an attracting fixed point for N (x). It turns out that the same is true if we replace dF (x) −1 by a fixed matrix C, which is sufficiently close to dF (x * ) −1 . The modified Newton operator is given by Now let us turn the things around and ask how can we use N m as a way to prove the existence of solution of (55). This is quite obvious. Namely, if U is homeomorphic to a closed finitedimensional ball and if then from the Brouwer Theorem it follows, that there exists To obtain the uniqueness it is enough show that N m is a contraction on U . Observe that it is impossible to verify in a single interval evaluation of the formula (57), that for some interval set ). It turns out the mean value form of N m can cure this deficiency.
This explains why the requirement has something to do with zeros of F (x).

The Krawczyk method
A method proposed by Krawczyk for finding zero's of F : . Typically x 0 is chosen to be midpoint of [x], we will denote this by ).
• C ∈ R n×n be a linear isomorphism The Krawczyk operator is given by Observe that point 2. in the above theorem gives us the way to establish the existence of unique zero of F in [x], while point 3. rules out the existence of zero in [x] i.e. in the terminology of previous section this is the exclusion test. When [x] is close to a zero of F then < F ([x]) > the evaluation of F ([x]) in the interval arithmetic might produce such overestimates that 0 ∈< F ([x]) >, while the Krawczyk operator will rule out the existence of 0 of F in [x]. This is in fact quite common phenomenon.
The Krawczyk operator is used as a part of iteration process In our context the only weakness of the Krawczyk operator requires the sets of the diameter smaller than 10 −2 to give us something. Above that threshold we usually have , F ) and the Krawczyk method is useless. This means that we still suffer with the curse of dimensionality.

The algorithm
The algorithm runs in the reduced configuration space which is a subset of R 2(n−1)−1 , i.e. a configuration is represented by a point (x 0 , y 0 , . . . , x n−3 , y n−3 , x n−2 ). Physically, we interpret such a configuration as n − 1 bodies with q i = (x i , y i ) for i = 0, . . . , n − 3 and q n−2 = (x n−2 , 0). From (44) we obtain q n−1 the position of the last body. Input: The input of the algorithm consists of 1. n -the number of bodies 2. some cube in the reduced configuration space.
Output: All different (up to reflections and rotations) central configurations in the full system for a given input cube. Since we use interval arithmetic for calculations, central configurations are also cubes containing the exact solution in their interior.
Program is divided into two stages: searching finds solutions and testing tests them to distinguish different CC and to find the kind of symmetry (if any exists).

Searching stage
In this stage we cover the configuration space with cubes. To fulfill requirements of Krawczyk's method (see Theorem 15) we must ensure that every point is in the interior of some cube. The algorithm runs for any initial cube, however if our goal is to find all the central configurations (for fixed n and equal masses) the reasonable cube is as follows (compare Sec. 6.1): Simple recursive algorithm works as follows: (I) if there is no solution in the cube return 0; (II) if there is unique solution in the cube return 1; (III) otherwise bisect the longest edge and recursively run the procedure for both parts; In the more detailed version of the algorithm the cube is represented by a vector of bodies. The class Bodies contains this vector and some methods to manipulate these bodies. An instance of Bodies is a cube bodies with some additional information.  (47) -(53)) 5. checkZero(bodies, i) -computes functions of vector field (see equations (4)) and tests if it is possible to have F (q 0 , . . . , q N −1 ) = 0.
To break a dimensionality curse we are looking for the possibility of restricting bodies before bisecting them (line 16). First place we are able to do this is the function thereIsNoSolution(bodies) (see Section 4.1). Another place is in krawczykMethod(bodies). If Krawczyk's method fails (line 5), bodies will have been restricted by intersection with the operator. Hence it is now (from line 17 onwards) the restricted cube that is being processed. This gives a large growth of efficiency.
Since the Krawczyk's method is costly and usually fails for large sets, we introduced a parameter bias. If the size (diameter) of all variables is not greater than bias then the Krawczyk's method is run. There is a big difference in execution time of the program depending on the value of the bias parameter; in Table 3 we present computation times for 5 bodies . For the same initial data program finds 8 solutions (some of them are later identified to be the same solutions), but numbers of failed and 'no-zero-inside' cubes differs.

Testing stage
The main goal of this stage is to identify distinct solutions. Additionally, we check the symmetry of solutions. In this stage we consider solutions in the full system. The solutions obtained in the search stage are given as a list boxes in which we have a unique solution. Some of these boxes may overlap and can in fact contain the same solution. Because we consider equal mass case we also do not want to distinguish solutions which differ by the indexes of the bodies. Hence two solutions produced in the searching stage can in fact be equivalent for two reasons: (1) the only difference is the ordering of the bodies, (2) the boxes defining them have non-empty intersection, having been obtained in different series of partitions. For any pair of solutions we treat the first one as a 'model', whilst bodies in the second solution are permuted in an attempt to match the model. When trying to match the solutions, we create a set containing both of them (an interval hull) and we prove the uniqueness within it. The rough idea is shown below: return (cZeros == 1); 10 } 11 } It may happen that there exists a solution in the set unionBodies but the set is too small to prove this using the Krawczyk operator, thus we 'inflate' it and retry the proof in the bigger set (the function blowUp(unionBodies)).
Establishing a (reflectional) symmetry of CC is very similar to testing the uniqueness of solutions. However, it requires an additional step to calculate a symmetric image. We take an interval hull of CC and its symmetric image. The function is as follows: 1 bool findSymByLine(Bodies &bodies) { 2 calculateParametersDefiningTheLine(p,q,a,b); //the line is (p,q) + t(a,b) 3 if (!isSingular(p * a + q * b)) //the line does not cross 0 4 return false; 5 Bodies symBodies, unionBodies; 6 calculateSymmetricImageOfTheSolution(bodies, p, q, a, b, symBodies)); 7 intervalHullBodies(bodies, symBodies, unionBodies); 8 int cZeros = search(unionBodies); 9 if (cZeros == 1) 10 return true; //bodies and symBodies are symmetric 11 else { 12 cZeros = blowUp(unionBodies); 13 return (cZeros == 1); 14 } First, in line 2, we calculate the parameters of the possible reflection symmetry line, but the symmetry tested contains also a permutation of bodies, we construct a configuration symBodies considering all possible permutations of bodies. Since the center of mass is in (0, 0), the line of reflection symmetry (p, q) + t(a, b) has to cross it, thus in line 3 we test this condition.
Note that lines 7-14 in the function findSymByLine are identical, up to the variable names, to lines 3-11 in theSameSolutions.

Technical data
The main computations were carried out in parallel using the template function std::async (from the standard C++ library) which runs the function asynchronously (potentially in a separate thread which may be part of a thread pool) on Dell R930 4x Intel Xeon E7-8867 v3 (2,5GHz, 45MB), 1024 GB RAM. The compiler is gcc version 4.9.2 (Debian 4.9.2-10+deb8u2). The best obtained times for different number of bodies with bias = 10 −2 are presented in Table 4  8 Minimizing dependency problem in gravitational force evaluation In this section we describe a method of calculations of functions and their derivatives used in the program. Let us denote and their derivatives are (with analogs for y's): In the program we have to evaluate f , where we just plug-in the interval arguments might lead to severe overestimation due to the dependency problem [Mo66, N90]. The best solution would be a cheap but rigorous estimate of sup and inf of f 1 i and f 2 i over D; this however appears to be a difficult and costly task. For the Krawczyk method we also need good estimates for ∂f i ∂x j and ∂f i ∂y j and we face the same problem. Thus we decided to optimize the computation of , since such components appear in all above equations. In equations (60) and (61) we evaluate expressions in the form xy/r 5 as a product of x/r 3 and y/r 2 which are treated as separated expressions.

8.1
Estimates for x a /r b and y a /r b Assume we want to calculate where We want to estimate f x and f y on D. We always assume that (0, 0) ∈ D. To minimize the overestimation of these calculations we look for the possible local extrema in D. Let us consider the function f x ; the second case of f y is symmetrical. First notice that by solving the system of equations we obtain (x, y) = (0, 0), which is impossible in our settings. Since there is no local extremum inside D thus it is attained on the edge of D. Te determine the specific points, where this extremum can be, we explore the equations (64) and (65) and analogical for f y , and we obtain: • border points on the lines y = ±x b−a a and x = ±y b−a a (blue points in Figure 4) • border points for x = 0 and for y = 0 (black points in Figure 4) Additionally we consider corners of D (red points in Figure 4). Next we examine all these point to establish the maximum and the minimum value of f x and f y . Note, that there is still a lot of room for further optimization, but for now only this version is implemented in the program. A Results and conclusions from the run of the program Below we give an output from our program (report files) for various number of bodies. These are proofs of Theorems 1 and 2. In the computations the masses were normalized so that m i = 1 and set λ = 1 in (2). Other authors have used different normalizations, both in [Moe89] and [F02] all masses were equal to 1, and in [Moe89] the condition I = 1 was demanded, while in [F02] the following normalization has been used |q 1 | |q 2 | · · · |q n−1 | = 1 |q n |, with q n−1 = (1, 0).
It is easy to see that the following quantity is an invariant for the transformations of space and mass scaling. We think that this is a meaningful invariant which makes comparison between CC obtained in different works, as it gives an obvious formula for rescaling of CC normalized in various ways.
In order to make a comparison with [Moe89, F02] easier we also include the Moeckel's potential, which equal to U (q) where q is CC scaled so that it satisfies I = 1 and m i = 1. This Moeckel's potential coincides with U √ I contained in With each CC we tried to prove that • it has a reflectional symmetry with respect OX-axis • it has a reflectional symmetry with respect to some other line. This other line is obtained as the perpendicular bisector of the segment joining q n−2 (the body furthest from 0) with other bodies (different bodies are tested).
Each of these reflections act in R 2 and are accompanied with a suitable permutation σ of bodies, which is displayed as a list of σ(0), . . . , σ(n − 1). For example, for n = 3 the list 2, 1, 0 means that σ(0) = 2, σ(1) = 1 and σ(2) = 0. The meaning of parameters in report files: • eps = ε, if the diameter of the box in max norm is smaller than this number then program stops subdividing this box and returns the message that it can not conclude whether there is CC in this box. This never happened in our computations • bias -the Krawczyk method is applied to a box only if the diameter of the box is smaller than bias.

A.1 Listing of central configurations
In this and nest sections we present all CCs for n = {3, . . . , 7} and asymmetrical CCs for n = {8, 9, 10} listed in [F02]. For any CC we give the picture of it (with red lines showing symmetry lines), an invariant J (see (66)), its coordinates and Moeckel's potential P . In the case of 3 and 4-body problem all CCs are symmetrical with respect to the 0X axis (this is the axis passing through the origin and the body most distant from the origin). In the case of 5, 6 and 7-body problems all CCs are still symmetrical with respect to some line, but there exists CCs which are not symmetrical with respect to the 0X axis. The first asymmetrical CCs appear for n = 8. Note that in the case of 6-body problem we have two different CCs which we called 'two equilaterals' and two which we called 'two isosceles triangles'. In Figures A.1.4.3 and A.1.4.7 a big triangle contains a small one whilst in Figures A.1.4.6 and A.1.4.8 there are two overlaping triangles. CCs in Figures A.1.4.5 and A.1.4.6 look almost the same, but in fact the first one is a real hexagon whist the second is very close to it, but it has not symmetry line in 0y -see Figure 5. In 7-body problem there two CCs, which look almost the same: the first one is a 'propeller'; the second is very close to it, but has only one symmetry line in 0x.   Moeckel's potential = [26.56875757, 26.56875757] permutation: 5, 2, 1, 3, 4, 0, symmetric with respect to OX no 6 permutation: 0, 4, 5, 3, 1, 2, reflectional symmetry with respect to other line 9 Number of different cc = 9 Seven bodies       2, 3, 4, 5, 6, 7, 8, 9, there is no symmetry with respect to OX Moeckel's potential = [118.3555704, 118.3555704] permutation: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, there is no symmetry with respect to OX