Stability Investigation of the Exact Symmetrical Solutions 
of the Plane Newton's Many-Body Problem

Alexander N. Prokopenya
Brest Technical University (Republic of Belarus)
E-mail: prokopenya@belpak.brest.by

Initialization

Default Settings

Clear["Global`*"];
SetOptions["stdout", 
        PageWidth -> 80];
$DefaultFont = {"Helvetica-Bold", 
          12};
th = AbsoluteThickness[2];

Off[General::spell1];
Off[General::spell];

Loading of Packages

Needs["Graphics`Colors`"];
Needs["Graphics`FilledPlot`"];

Introduction

The many-body problem was first formulated precisely by Newton. In the simplest case when the objects involved are point masses it may be stated as follows: Given at any time the positions and velocities of three or more massive particles moving under their mutual gravitational forces, the masses also being known, calculate their positions and velocities for any other time. The problem is more complicated when the bodies' shapes and internal constitutions have to be taken into account as in the Earth-Moon-Sun problem.Although this problem is very important for the wide variety of applications, only in the case of two interacting particles the corresponding general solution is found. Moreover, it turned out the differential equations of the many-body problem are in general not integrable. So according to Poincare ideas the further progress in this field will be connected with finding new classes of the exact particular solutions of the many-body problem and investigating their stability. The stability investigation turned out to be the most complicated problem of qualitative theory of differential equations because it requires a lot of tremendous analytical calculations. For example, solving of the stability problem of the Lagrange triangular solutions has taken about 200 years, whereas stability of the homographic and homothetic solutions in the three-body problem still remains unsolved. However, now we have new computation systems, for example Mathematica, that essentially increases our ability in doing both numerical and symbolic calculations. So there is a hope that we can push considerably both the many-body problem and the theory of dynamical systems generally using modern computer algebra systems.
	In [1] B. Elmabsout found a new exact particular solution of the plane Newton’s many-body problem where the n bodies  (Subscript[P, 1]), ... , (Subscript[P, n]) of equal masses are at the vertices of a regular polygon with n sides that rotates uniformly around the body (Subscript[P, 0])of a non-zero mass (Subscript[m, 0]) being at the center of polygon. In the case of n=5 the corresponding configuration is shown on Fig. 1.

Then E.A.Grebenikov showed [2, 3] that there exist a whole class of the exact symmetrical solutions of the problem of (n+1)-bodies that are represented by regular, scale-similarly to self varying polygons, rotating with variable angular velocity around their centers (Fig. 2).

Stability of the invariant polygon was completely investigated in linear approximation in [4]. So the problem of studying stability of the exact symmetrical solutions is very topical and interesting.

I. Equations of Motion in Nechvil's variables

Let us consider the (n+1) particles (Subscript[P, 0]), (Subscript[P, 1]), ... , (Subscript[P, n]) of masses (Subscript[m, 0]), (Subscript[m, 1]) = (Subscript[m, 2]) =  ... (Subscript[m, n]) = 
            m moving under their mutual gravitational forces (Subscript[F, ij]) = 
        f (((Subscript[m, i]) (Subscript[m, j]))/((Subscript[Δ, ij])^2)), where (Subscript[Δ, ij]) is the distance between i-th and j-th particles and f is a gravity constant. In a relative Cartesian frame (Subscript[P, 0])xyz with the origin being at (Subscript[P, 0]) the motion equations of the particles (Subscript[P, 1]), ... , (Subscript[P, n]) are

\(TextForm\`\(\(d\^2\) \(r\& -> \)\_j\)\/\(d\ t\^2\) + 
        f \((m\_0 + m)\)\ \(r\& -> \)\_j\/r\_j\^3 = 
      f\ m \(∑\+\(i = 1\ \((\(≠\)\(j\))\)\)\%n\((\(\(r\& -> \)\_i - \(r\& -> \)\_j\)\/\(\(|\)\(\(r\& -> \)\_i - \(r\& -> \)\_j\)\( | \^3\)\) - \(r\& -> \)\_i\/r\_i\^3)\)\ \ \ \ \ \ \ \ \ \ \ \ \(\((j = \(1, n\)\&_)\)\(.\)\)\)\)

Evaluating the sum in the right part of (1) we should exclude the case of i = j, so we can not use the built-in function Sum. But we can replace the sum on i with a multiplier Sumi and handle equations (1) without any problems. Of course we should be careful doing some transformation of equations (1) to avoid the situation when Sumi turns out to be in the denominator.  Thus, we define equations (1) as eq1 using the next notations:
r[i] = (Sqrt[((Subscript[x, i])^2) + ((Subscript[y, i])^2) + ((Subscript[z, i])^2)]) - a distance between the particle (Subscript[P, i]) and the origin;
\(TraditionalForm\`R[i, 
          j] = \@\(\((x\_i - x\_j)\)\^2 + \((y\_i - y\_j)\)\^2 + \((z\_i - z\_j)\)\^2\)\)  - a distance between two particles (Subscript[P, i]) and (Subscript[P, j]), for i,j=\(\(1, n\)\&_\).
Projections of the expressions under the sum in (1) we denote as Xj[i, j], Yj[i, j], Zj[i, j] and define as a rule  r1.

r1 = {Xj[i, j] -> (((Subscript[x, i]) - (Subscript[x, j]))/R)[i, j]^3 - ((Subscript[x, i])/r)[i]^3, 
          Yj[i, j] -> (((Subscript[y, i]) - (Subscript[y, j]))/R)[i, j]^3 - 
              ((Subscript[y, i])/r)[i]^3,  
    Zj[i, 
              j] -> (((Subscript[z, i]) - (Subscript[z, j]))/R)[i, j]^3 - 
              ((Subscript[z, i])/r)[i]^3};

eq1 = (∂_(t, t)){(Subscript[x, j])[t], (Subscript[y, j])[t], 
              (Subscript[z, j])[t]} == ({-f ((Subscript[m, 0]) + m) ((Subscript[x, j])/r)[j]^3 + 
                    f m Sumi Xj[i, 
                        j], 
                                                         -f ((Subscript[m, 0]) + m) ((Subscript[y, j])/r)[j]^3 + 
                    f m Sumi Yj[i, 
                        j], 
                                                          -f ((Subscript[m, 0]) + m) ((Subscript[z, j])/r)[j]^3 + 
                    f m Sumi Zj[i, 
                        j]} /. 
                                               r1  /.  {(Subscript[x, i])_ -> (Subscript[x, i])[t], (Subscript[y, i])_ -> (Subscript[y, i])[t], 
                (Subscript[z, i])_ -> 
                  (Subscript[z, i])[t]})  // 
                                                                        Thread[#, List] &

It is convenient to rewrite the equations of motion eq1 using cylindrical coordinates in accordance with the relationship \(TraditionalForm\`\(r\& -> \)\_j = \((ρ\_j\ cosφ\_j, \ ρ\_j\ sinφ\_j, \ z\_j)\)\). Defining a transformation rule r2 we obtain the equations eq2.

r2 = {(Subscript[x, j])_ -> ((Subscript[ρ, j])[#] Cos[(Subscript[φ, j])[#]] &), 
          (Subscript[y, j])_ -> ((Subscript[ρ, j])[#] Sin[(Subscript[φ, j])[#]] &), 
          (Subscript[z, j])_ -> ((Subscript[z, j])[#] &)} ;
eq2 = eq1  /. r2 // Simplify;

Solving the first and the second equations of the system eq2 with respect to (Subscript[ρ, j])'' (t) and (Subscript[φ, j])'' (t) we rewrite eq2 as the system eq3.

sol2 = Solve[{eq2[[1]], 
                eq2[[2]]}, {(Subscript[ρ, j])''[t], (Subscript[φ, j])''[
                  t]}][[1]]  // 
          Simplify;
eq3 = {(Subscript[ρ, j])''[
            t] == ((Subscript[ρ, j])''[t] /. sol2  // 
              Collect[#, Sumi, 
                  Simplify] &), 
        (Subscript[ρ, j])[t] (Subscript[φ, j])''[
              t] == ((Subscript[ρ, j])[t] (Subscript[φ, j])''[t] /. sol2  // 
              Simplify), 
          eq2[[3]] }

A distance r[i] between the particle (Subscript[P, i]) and the origin in cylindrical coordinates is defined in a rule r3.

r3 = 
      r[j_] -> ((Sqrt[(Subscript[x, j])[t]^2 + (Subscript[y, j])[t]^2 + (Subscript[z, j])[t]^2])  /. r2  // 
            Simplify)

A distance between two particles (Subscript[P, i]) and (Subscript[P, j]) is defined in a rule r4.

r4 = 
      R[i_, j_] -> ((Sqrt[((((Subscript[x, i])[t] - (Subscript[x, j])[t]))^2) + ((((Subscript[y, i])[t] - (Subscript[y, j])[t]))^2) + ((((Subscript[z, i])[t] - (Subscript[z, j])[t]))^2)])  /. r2  // 
            Simplify)

It is easy to verify that equations eq3 have a solution

\(TextForm\`\(ρ\_j\) \((t)\) = \(ρ\_1\) \((t)\) ≡ ρ \((t)\), \ \(φ\_j\) \((t)\) = φ \((t)\) + \(\(2  π\)\/n\) j, \ \(z\_j\) \((t)\) = 
      0\ \ \ \ \ \ \ \ \((j = \(1, n\)\&_)\)\)

where the functions ρ(t) and φ(t) are connected with two relationships

\(TextForm\`ρ \((φ)\) = 
      p\/\(1 + e\ cosφ\), \ \ \ ρ\^\(\(\ \)\(2\)\)\ \(d\ φ\)\/\(d\ t\) = \(c = \(\(const\)\(\ \)\(.\)\)\)\)

Parameter p is determined with the constant c according to the relationship

p = ((c^2)/(f (Subscript[M, n]))), where  (Subscript[M, n])=(Subscript[m, 0])+(m/4)(∑_(j = 1)^(n - 1))(1/(Sin ((π/n) j))).

So two parameters p and e completely determine trajectories of the particles. In the case of e = 0 the particles move uniformly on the circle of radius p  with angular velocity ω that equals

ω = (Sqrt[((f (Subscript[M, n]))/(p^3))]).

To investigate stability of the solution (2) it is convenient to transform equations eq3 to Nechvil's phase space according to the rule r5.

r5 = {(Subscript[ρ, j])_[
              t] :> ((p/(1 + e Cos[φ])) (Subscript[ρ, j])[φ]), 
          (Subscript[z, j])_[t] :> ((p/(1 + e Cos[φ])) (Subscript[z, j])[φ]), (Subscript[φ, j])_[
              t] :> (Subscript[φ, j])[φ]};

where the polar angle φ is considered as a new independent variable. Using the relationship

\(TextForm\`ρ\^2\ \(d\ φ\)\/\(d\ t\) = \(p\^2\/\((1 + e\ cosφ)\)\^2\ \(d\ φ\)\/\(d\ t\) = \(c = \(\@\(p\ f\ M\_n\) = const\)\)\)\)

we can define a rule for derivatives transformation as

r6 = {D[f_[t], {t, n_}] :> 
            Nest[(c/(p^2)) (((1 + e Cos[φ]))^2) D[#, φ] &, f[t] /. r5, n ]};

Applying the rule r6 to derivatives of the functions (Subscript[ρ, j])(t), (Subscript[φ, j])(t), (Subscript[z, j])(t) we obtain the next results.

(Subscript[ρ, j])'[t] /. r6  // Simplify

(Subscript[ρ, j])''[t] /. r6  // Simplify

(Subscript[φ, j])'[t] /. r6  // Simplify

(Subscript[φ, j])''[t] /. r6  // Simplify

(Subscript[z, j])'[t] /. r6  // Simplify

(Subscript[z, j])''[t]  /. r6  // Simplify

Nechvil's transformation means that all distances in equations of motion get a multiplier (p/(1 + e cosφ)). So the distances r[j] and R[i,j] in eq3 should be transformed according to the rules r7 and r8.

r7 = r[j_] -> (p/(1 + e Cos[φ])) r[j];

r8 = 
        R[i_, j_] -> 
          (p/(1 + e Cos[φ])) R[i, j];

Taking into account the rules r5, r6, r7 and r8 we can rewrite equations eq3 as

eq4 = 
      Map[Collect[#, Sumi, 
            Simplify] &, {(Subscript[ρ, j])''[φ] == ((Subscript[ρ, j])''[φ]  /. Solve[eq3[[1]]  /. r5  /. r6  /. r7  /. 
                      r8, (Subscript[ρ, j])''[φ]][[1]]) , 
(Subscript[ρ, j])[φ] (Subscript[φ, j])''[φ] == ((Subscript[ρ, j])[φ] (Subscript[φ, j])''[φ] /. Solve[eq3[[2]]  /. r5  /. r6  /. r7  /. 
                      r8, (Subscript[φ, j])''[φ]][[1]]), 
(Subscript[z, j])''[φ] == ((Subscript[z, j])''[φ]  /. Solve[eq3[[3]]  /. r5  /. r6  /. r7  /. 
                      r8, (Subscript[z, j])''[φ]][[1]])}, 2]

In Nechvil's variables solution (2) can be written as

\(TextForm\`\(ρ\_j\) \((φ)\) = ρ ≡ 
        1, \ \(φ\_j\) \((φ)\) = φ + \(\(2  π\)\/n\) j, \ \(z\_j\) \((φ)\) = 
      0\ \ \ \ \ \ \ \ \((j = \(1, n\)\&_)\)\)

To verify whether functions (4) are the solution of the equations of motion let us substitute them into equations eq4. The corresponding substitutions are defined in rules r9, r10, r11.

r9 = 
        r[j_] -> (r[j] /. r3 /. 
              t -> φ);
r10 = 
        R[i_, j_] -> (R[i, j] /. r4 /. 
              t -> φ);
r11 = {(Subscript[ρ, j])_ -> (1 &), (Subscript[φ, j])_ -> (# + ((2  π)/n) j &), (Subscript[z, j])_ -> (0  &)};

eq5 = eq4  /. r9  /. r10  /. r11  // Simplify

There are the next sums in equations eq5:

(∑_(i = 1 (≠j))^n) 
        Cos[((2  π)/n) (i - j)],

(∑_(i = 1 (≠j))^n) 
        Sin[((2  π)/n) (i - j)],

\(\(\(∑\+\(i = 1 \((\(≠\)\(j\))\)\)\%n\(\(\ \)\(Cot[\(π\/n\) \((i - j)\)]\)\)\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\(|\)\)\)\(,\)\)\)

(∑_(i = 1 (≠j))^n) 
            (1/(|Sin[(π/n) (i - 
                      j)]|)).

It is easy to prove that these sums do not depend on j. For example, the first sum may be rewritten as

\(∑\+\(i = 1 \((\(≠\)\(j\))\)\)\%n 
          Cos[\(\(2  π\)\/n\) \((i - 
                j)\)] = 
\(\(=\)\(∑\+\(i = 1\)\%\(j - 1\)Cos[\(\(2  π\)\/n\) \((i - j)\)] + ∑\+\(i = j + 1\)\%n 
                Cos[\(\(2  π\)\/n\) \((i - 
                      j)\)] = 
\(\(=\)\(∑\+\(i = 1\)\%\(j - 1\)Cos[\(\(2  π\)\/n\) \((n + i - 
                          j)\)] + ∑\+\(i = 1\)\%\(n - j\)Cos[\(\(2  π\)\/n\) 
                      i] = 
\(\(=\)\(∑\+\(i = n - j + 1\)\%\(n - 1\)Cos[\(\(2  π\)\/n\) 
                          i] + ∑\+\(i = 1\)\%\(n - j\)Cos[\(\(2  π\)\/n\) i] = 
\(\(=\)\(∑\+\(i = 1\)\%\(n - 1\)\(\(Cos[\(\(2  π\)\/n\) i]\)\(.\)\)\)\)\)\)\)\)\)\)\)

The sum obtained may be calculated easily.

Sum[Cos[((2 π)/n) k], {k, 1, n - 1}] // FullSimplify

We can easily verify that even without any transformations the sum may be calculated. Actually, as in the case of i=j the corresponding term Cos[(2/n) (i - j)] = 1 then subtracting this term in the sum we do not change the result. So we have

Sum[Cos[((2  π)/n) (i - j)],  {i, n}] - 
      Cos[((2  π)/n) (j - j)]

We can calculate the second some in the similar way.

Sum[Sin[((2  π)/n) (i - j)],  {i, n}]

The third sum can not be calculated with Mathematica but it may be easily shown to equal zero.

\(\(∑\+\(i = 1 \((\(≠\)\(j\))\)\)\%n Cot[\(π\/n\) \((i - j)\)]\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\(|\)\) = \(∑\+\(i = 1 \((\(≠\)\(j\))\)\)\%n\( Cos[\(π\/n\) \((i - j)\)]\ Sin[\(π\/n\) \((i - j)\)]\)\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\( | \^3\)\) = 
\(\(=\)\(∑\+\(i = 1\)\%\(j - 1\)\(Cos[\(π\/n\) \((i - j)\)]\ Sin[\(π\/n\) \((i - j)\)]\)\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\( | \^3\)\) + ∑\+\(i = j + 1\)\%n\( Cos[\(π\/n\) \((i - j)\)]\ Sin[\(π\/n\) \((i - j)\)]\)\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\( | \^3\)\) = 
\(\(=\)\(∑\+\(i = 1\)\%\(j - 1\)\(Cos[\(π\/n\) \((n + i - j)\)]\ Sin[\(π\/n\) \((n + i - j)\)]\)\/\(\(|\)\(Sin[\(π\/n\) \((n + i - j)\)]\)\( | \^3\)\) + ∑\+\(i = 1\)\%\(n - j\)\(Cos[\(π\/n\) i]\ Sin[\(π\/n\) i]\)\/\(\(|\)\(Sin[\(π\/n\) i]\)\( | \^3\)\) = 
\(\(=\)\(∑\+\(i = n - j + 1\)\%\(n - 1\)\(Cos[\(π\/n\) i]\ Sin[\(π\/n\) i]\)\/\(\(|\)\(Sin[\(π\/n\) i]\)\( | \^3\)\) + ∑\+\(i = 1\)\%\(n - j\)\(Cos[\(π\/n\) i]\ Sin[\(π\/n\) i]\)\/\(\(|\)\(Sin[\(π\/n\) i]\)\( | \^3\)\) = 
\(\(=\)\(∑\+\(i = 1\)\%\(n - 1\)\(Cos[\(π\/n\) i]\ Sin[\(π\/n\) i]\)\/Sin[\(π\/n\) i]\^3 = \(∑\+\(i = 1\)\%\(n - 1\)Cos[\(π\/n\) i]\/Sin[\(π\/n\) i]\^2 = 
                            0\)\)\)\)\)\)\)\)\)\);\)\)

The fourth sum is a positive number that we denote (Subscript[S, 1]).

\(∑\+\(i = 1 \((\(≠\)\(j\))\)\)\%n 1\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\(|\)\) = \(∑\+\(i = 1\)\%\(j - 1\)1\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\(|\)\) + ∑\+\(i = j + 1\)\%n 1\/\(\(|\)\(Sin[\(π\/n\) \((i - j)\)]\)\(|\)\) = 
\(\(=\)\(∑\+\(i = n - j + 1\)\%\(n - 1\)1\/\(\(|\)\(Sin[\(π\/n\) i]\)\(|\)\) + ∑\+\(i = 1\)\%\(n - j\)1\/\(\(|\)\(Sin[\(π\/n\) i]\)\(|\)\) = 
\(\(=\)\(∑\+\(i = 1\)\%\(n - 1\)1\/Sin[\(π\/n\) i] ≡ \(\(S\_1\)\(.\)\)\)\)\)\)\)\)

Taking into account the results obtained we see that the second equation in eq5 is true. Substituting the corresponding sums into the first equation of eq5 we obtain the equation eq6.

eq6 = Map[Expand, eq5[[1]], 
              1]  /. Sumi Cos[((2 (i - j) π)/n)] -> -1  /. 
          Sumi Csc[(((i - j) π)/n)]^2 (Sqrt[Sin])[(((i - j) π)/n)]^2 -> (Subscript[S, 1])  // Simplify

It is easy to see that the equation obtained is true if accorging to (3) we use the equality (c^2) = f p (Subscript[M, n]) and (Subscript[M, n]) = ((Subscript[m, 0]) + (m/4) (Subscript[S, 1])).

eq6   /. (c^2) ->  p f (Subscript[M, n]) /. 
        (Subscript[M, n]) -> ( (Subscript[m, 0]) + (m/4) (Subscript[S, 1]))  // Simplify

So we see that the functions (2) and (4) are the solutions of the equations of motion that are determined in Nechvil's variables as eq4.

II. Linearization of the Equations of Motion

We shall study the equations eq4 at the neighbourhood of the solution (4). So we should make the next replacement in eq4:

\(TextForm\`ρ\_j -> 
      1 + u\_j, \ φ\_j -> φ + \(2  π\)\/n\ j + γ\_j, \ j = \(1, n\)\&_\)

considering the functions (Subscript[u, j])(φ), (Subscript[γ, j])(φ), (Subscript[z, j])(φ) as small perturbations of the solutions (4). To obtain linearized equations of motion it is necessary to expand equations eq4 in powers of (Subscript[u, j]), (Subscript[γ, j]) and (Subscript[z, j]) neglecting all terms of order superior or equal to 2. To do this in a simplest way we supply the functions (Subscript[u, j])(φ), (Subscript[γ, j])(φ), (Subscript[z, j])(φ) with a multipier δ according to the rule r12.

r12 = {(Subscript[ρ, j])_ -> (1 + δ (Subscript[u, j])[#] &), (Subscript[φ, j])_ -> (# + ((2  π)/n) 
                    j + δ (Subscript[γ, j])[#] &),  
          (Subscript[z, j])_ -> (δ (Subscript[z, j])[#] &)};

Then equations eq4 may be rewritten as eq7.

eq7 = eq4  /. r9  /. r10  /. r12;

Now it is enough to expand equations eq7 in powers of δ neglecting all terms of order superior 1. As a result we obtain the next system of differential equations:

linearization[
          f_Equal] := Series[f, {δ, 0, 1}]  // 
            Thread[#, Equal] &  // Normal ;
eq8 = linearization[#] &  /@  eq7;

Of course in the case of  δ -> 0 equation eq8 should coincide with equations eq5 obtained above. Actually, making this replacement we obtain the equations eq5.

eq8 /. δ -> 0  // Simplify

Extracting coefficients of δ in each equation of eq8 we obtain a linearized system of differential equations.

linearApproximation[f_Equal] := 
        Coefficient[f[[1]], δ] - 
            Coefficient[
              f[[2]], δ] == 
          0;
eq9 = linearApproximation[#] &  /@  eq8;

Solving the system eq9 with respect to (Subscript[u, j])(φ), (Subscript[λ, j])(φ) and (Subscript[z, j])(φ) we can rewrite linearized equations of motion as eq10, eq11 and eq12. Using the next notations

(c^2) = f p (Subscript[M, n]) = f p ((Subscript[m, 0]) + (m/4) (Subscript[S, 1])) = 
          f m p (μ + (1/4) (Subscript[S, 1])),  μ = 
      ((Subscript[m, 0])/m)

we obtain

eq10 = (Subscript[u, j])''[φ]==((Subscript[u, j])''[φ] /. Solve[eq9[[1]]  /. 
                      c -> (Sqrt[f m p (μ + (1/4) (Subscript[S, 1]))]) /. 
                    (Subscript[m, 0]) ->  
                      m μ, (Subscript[u, j])''[φ]][[1]]  // 
            Collect[#, Sumi, Simplify] &)

eq11 = (Subscript[γ, j])''[φ]==((Subscript[γ, j])''[φ] /. Solve[eq9[[2]]  /. c -> (Sqrt[f m p (μ + (1/4) (Subscript[S, 1]))]) /. 
                    (Subscript[m, 0]) ->  
                      m μ, (Subscript[γ, j])''[φ]][[1]]  // 
            Collect[#, Sumi, Simplify] &)

eq12 = (Subscript[z, j])''[φ]==((Subscript[z, j])''[φ] /. Solve[eq9[[3]]  /. 
                      c -> (Sqrt[f m p (μ + (1/4) (Subscript[S, 1]))]) /. 
                    (Subscript[m, 0]) ->  
                      m μ, (Subscript[z, j])''[φ]][[1]]  // 
            Collect[#, Sumi, Simplify] &)

Now the stability problem of solution (4) in linear approximation is reduced to studying stability of a trivial solution of the system eq10,eq11, eq12 that is a linear system of differential equations with periodic coefficients. So we should find and analyze the characteristic exponents of this system.

III. Studying Stability of the Solutions for Perpendicular Perturbations

From the system above we can see that functions (Subscript[z, j]) determining perturbations of the particles' orbits in perpendicular to the plane (Subscript[P, 0])xy direction do not depend on (Subscript[u, j]) and (Subscript[γ, j]). So in linear approximation we can study stability of the solutions (4)  in respect to the perpendicular perturbations (Subscript[z, j]) and the perturbations (Subscript[u, j]) and (Subscript[γ, j]) contained in the plane of orbits (Subscript[P, 0])xy  separately. 
	Let us start with studying linear stability of the solutions for perpendicular perturbations. It is convenient to rewrite the system eq12 in vector notations. Thus we introduce the n-dimensional vector Z with components (Subscript[z, j])(φ) so that

Z=(({{(Subscript[z, 1])}, {   ⋮

                  (Subscript[z, n])}})),   ((d Z)/(d φ))=(({{(Subscript[(z'), 1])}, {   ⋮

                    (Subscript[(z'), n])}})).

We also introduce n×n matrixes

\(B : \ \ B\_jk = 
        1\/\(\(|\)\(Sin \((\(π\/n\) \((k - j)\))\)\)\( | \^3\)\)\ \((j ≠ k)\)\), \(B\_jj = \(-\(∑\+\(k = 1 \((\(≠\)\(j\))\)\)\%n b\_jk\)\)\),
K :  (Subscript[K, jk]) = 1,
I :  (Subscript[I, jk])=(Subscript[δ, jk])={({{1,  j = k}, {0,  j ≠ k}})   \(\((j, 
                k = \(1, n\)\&_)\)\(.\)\)

Then we can rewrite equations eq12  as

\(TextForm\`\(\(\ \ \)\(\(d\^2\ Z\)\/\(d\ φ\^2\) = \(1\/\(2 \((1 + e\ cosφ)\) \((4  μ + S\_1)\)\)\) \((B - 8  K - 
            8  μ\ I - 
            2 \((4  μ + 
                  S\_1)\)\ e\ cosφ\ I)\)\ \(\(Z\)\(.\)\)\)\)\)

It's easy to verify that the matrix B has eigenvectors (Subscript[E, r]) with components

\(\(\(E\_\(k, r\) = 
      1\/\@n\ Exp[i\ \(2  π\ r\)\/n\ \ k]\ \ \ \((k, 
          r = \(1, n\)\&_)\)\)\(,\)\)\)

where i = (Sqrt[-1]) and the first index k denotes the number of the component. The corresponding eigenvalues are (-(Subscript[λ, r])) where

\(λ\_r = 
      2 \(∑\+\(k = 1\)\%\(n - 1\)\(\(\(\(Sin\^2\) \((\(π\/n\) 
                      r\ k)\)\)\/\(\(Sin\^3\) \((π\/n\ k)\)\)\)\(.\)\)\)\)

Vectors (Subscript[E, r]) are also the eigenvectors of the matrix K that has two eigenvalues: 0 and n. Actually,

K . (Subscript[E, r])=0 for r=\(\(1, n - 1\)\&_\)     and K . (Subscript[E, n])=n (Subscript[E, n])

So, introducing the matrix Q with elements

\(Q : \ \ Q\_jk = 
      1\/\@n\ Exp[
          i\ \(2  π\ k\)\/n\ \ j]\ \ \ \(\((j, 
            k = \(1, n\)\&_)\)\(.\)\)\)

we can diagonilize the matrixes B and K according to the transformation rules

B  -> (Q^( +)) 
        B Q,   K  -> (Q^( +)) K Q

So we have reduced equation (5) to the normal form

\(\(\(\(\(d\^2\)\(\ \)\(z\_r\)\(\ \)\)\/\(d\ φ\^2\) = \(-\(\(a\_r + e\ Cosφ\)\/\(1 + 
                e\ Cosφ\)\)\)\ z\_r\)\(,\)\)\)

where

(Subscript[a, r]) = (((Subscript[λ, r]) + 8  μ + 
              8  n (Subscript[δ, r, n]))/(2 (4  μ + 
                  (Subscript[S, 1])))).

Equation (6) is the second order differential equation that is sometimes called a Hill's equation. In the case of e = 0 it reduces to the equation of simple harmonic oscillator of frequency (Sqrt[a]).
Proposition.  As (Subscript[a, r]) > 0 for any \(TraditionalForm\`r = \(1, n\)\&_\) in the case of e = 0 the solutions (4) are linearly stable  for any perpendicular to the plane of the particles' orbits perturbations.

Studying Stability of a Zero Solution of Hill's Equation

Let us investigate stability of a zero solution of the next Hill's equation.

eq20 = Z''[φ] == -((a + 
                    e Cos[φ])/(1 + 
                    e Cos[φ])) Z[φ];

It is apparent that equation eq20 can be easily reduced to a linear system of two first order differential equations with periodic matrix:

((d Z)/(d φ)) = 
      A (φ) Z,

where

Z=(({{(Subscript[z, 1])}, {(Subscript[z, 2])}})),  ((d Z)/(d φ))=(({{((d (Subscript[z, 1]))/(d φ))}, {((d (Subscript[z, 2]))/(d φ))}})),  A (φ)=(({{0, 1}, {-((a + e Cos[φ])/(1 + 
                        e Cos[φ])), 0}}))

General analysis of such systems shows that depending on the values of parameters a and e its zero solution may be stable or unstable [5]. The boundaries between the domains of stability and instability in the  a-e plane are some curves a = a(e) which are characterized by the presence of periodic solutions with period 2π or 4π. It is apparent that if some function Z=Z(φ) is a solution of eq20 then the functions Z = Z(-φ) and Z = -Z(φ) are its solutions too. Thus, we can conclude that there are even and odd periodic solutions of eq20. So we shall seek the even periodic solutions of period 2π  in the form of a Fourier series

\(\(Z = 
        c\_0 + ∑\+\(k = 1\)\%∞ c\_\(1  k\)\ cos \((k\ φ)\)\)\(,\)\)

and the odd periodic solutions of the same period in the form

\(Z = ∑\+\(k = 1\)\%∞ c\_\(2  k\)\ sin \(\((k\ φ)\)\(.\)\)\)

Periodic solutions of period 4π may be represented as

\(Z = 
        c\_0 + ∑\+\(k = 1\)\%∞ c\_\(1  k\)\ cos \((k\/2\ φ)\), \ \ Z = ∑\+\(k = 1\)\%∞ c\_\(2  k\)\ sin \((k\/2\ φ)\)\)

It's apparent that the last formulas can also be used to obtain the solutions with period 2π by setting to zero the Fourier coefficients corresponding to k being an odd integer.

The Even Solutions.

It is convenient to rewrite equation eq20 as

eq20 = 
      eq20 (1 + e Cos[φ])  // Thread[#, Equal] &

Substituting the third Fourier series into eq20 we obtain

eq21 = 
      eq20 /. Z -> (c0 + 
                  Sumk (c1[2  k - 1] Cos[((2  k - 1)/2) #] + 
                        c1[2  k] Cos[k #]) &)  // Simplify

Here we separated the terms with even and odd k. We also replaced the sign of sum with the corresponding multiplier Sumk. Reducing trigonometrical functions  according to the identities

Cos[α] Cos[β]  // TrigReduce

Sin[α] Cos[β]  // TrigReduce

we have found

reducing[
          f_Equal] := (f[[1]] - 
                  f[[2]]  // 
                Expand  // TrigReduce  // 
            Expand);
eq22 = reducing[ eq21 ]

Since this Fourier series must vanish for all φ, every coefficient must vanish. First, the constant term in the series must be equal to zero.

eq23 = (eq22  /. {k -> 1, 
                Sumk -> 1}  /. {Cos[_] -> 0, 
              Sin[_] -> 0}) == 0

Coefficients of  Cos[φ] in eq22 must be equal to zero too.

eq24 = (Coefficient[eq22, Cos[φ]] + 
            Coefficient[Coefficient[eq22 /. k -> 1, Sumk], 
              Cos[φ]] + 
            Coefficient[Coefficient[eq22 /. k -> 2, Sumk], 
              Cos[φ]]) == 0

Extracting coefficients of  Cos[k φ] in eq22 we obtain an infinite sequence of equations for k >1.

coeffCos1[f_] := 
        Coefficient[f, 
            Cos[k φ]] + (Coefficient[f, 
                Cos[φ - k φ]]  /. 
              k -> k + 1) + (Coefficient[f, 
                Cos[φ + k φ]] /. 
              k -> k - 1) ;

eq25 = (coeffCos1[eq22]  /. Sumk -> 1  // 
            Collect[#, {c1[-1 + 2  k], c1[2  k], c1[1 + 2  k]}, 
                Simplify] &) == 0

Coefficients of the Fourier series c1[2  k] are determined as a solution of this infinite sequence of linear homogeneous equations. As there is a periodic solution of the equation eq20 the system of equations obtained must has non-trivial solution. So the corresponding determinant of the system must equal zero, thus giving us an equation determining the stability boundaries a = a(e) in the a-e plane. To calculate the determinant we first define a function mat1(n) determining n × n matrix of the system.

mat1[n_, eq0_Equal, eq1_Equal, eq2_Equal] := 
        Module[{j, l, 
            m}, 
                m = {Join[{Coefficient[eq0[[1]], c0]}, 
                Table[Coefficient[
                    eq0[[1]], 
                    c1[2  j]], {j, 
                    n}]]}; 
                 If[n > 0,  
            m = Append[m, 
                Join[{Coefficient[
                      eq1[[1]], c0]}, 
                  Table[Coefficient[
                      eq1[[1]], 
                      c1[2  j]], {j, 
                      n}]]]]; 
                 If[n > 1, 
            Do[m = Append[m, 
                  Join[{Coefficient[
                        eq2[[1]], c0]}, 
                    Table[Coefficient[
                        eq2[[1]] /. 
                          k -> l, c1[2  j]], {j, 
                        n}]]], 
                              {l, 2, 
                n}]]; 
                           m ];

In the case of n=6 for example the corresponding matrix is

mat1[6, eq23, eq24, eq25]  // MatrixForm

If we equalize determinant of this matrix to zero for n -> ∞ we'll obtain an equation giving us a function a = a(e) for which the even periodic solution of eq20 exist. Of course it's impossible to calculate a determinant of the infinite matrix. So let us calculate the determinant of the matrix for a finite n. It is easy to see that the next recurrence relation is true for the determinants.

d1[0] = a;
d1[1] = a (a - 1);
d1[k_] := (a - (k^2)) d1[k - 1] - 
          ((e^2)/4) (k - 2) (k - 1) k (k + 1) d1[
              k - 2];

From the formula above it's apparent that in the case of e = 0 determinants will equal to zero when a = (k^2)  (k = 0,  1,  2,  ...  ). It means the stability boundaries should cross the e = 0 axis in the a e-plane at these points. For sufficiently small e we can represent the stability boundaries in the vicinity of the points a = (k^2) as a power series in e. Let us calculate, for example, the determinant in the vicinity of the point a = (2^2) truncating the the corresponding power series after the term of the fourth order.

eq26 = 
      d1[3] /. a -> 4 + (Subscript[δ, 1]) (e^2) + (Subscript[δ, 2]) (e^4)  // 
        Collect[#, e, Simplify] &

Setting the determinant to zero we find the coefficients

eq27 = {Coefficient[eq26, e, 2] == 0,  
          Coefficient[eq26, e, 4] == 0};
sol27 = Solve[
          eq27, {(Subscript[δ, 1]), (Subscript[δ, 2])}][[1]]

In order to improve the results we increase the order of the determinant.

eq28 = 
        d1[4] /. a -> 4 + (Subscript[δ, 1]) (e^2) + (Subscript[δ, 2]) (e^4)  // 
          Collect[#, e, Simplify] &;

eq29 = {Coefficient[eq28, e, 2] == 0,  
          Coefficient[eq28, e, 4] == 0};
sol29 = Solve[
          eq29, {(Subscript[δ, 1]), (Subscript[δ, 2])}][[1]]

Taking into account determinant of the next order we obtain

eq30 = 
        d1[5] /. a -> 4 + (Subscript[δ, 1]) (e^2) + (Subscript[δ, 2]) (e^4)  // 
          Collect[#, e, Simplify] &;

eq31 = {Coefficient[eq30, e, 2] == 0,  
          Coefficient[eq30, e, 4] == 0};
sol31 = Solve[
          eq31, {(Subscript[δ, 1]), (Subscript[δ, 2])}][[1]]

It's apparent from the recurrence relation that if we want to find the stability boundary a = a(e) crossing the axis e = 0 in the point a = (k^2)with an accuracy (e^(2  p)) we should use determinant of the (k+p)-th order. As a result we obtain the curves

ae1[k_] := 
        Module[{q1, q2, δ1, δ2, s1, k1},  
          k1 = d1[k + 
                2]; 
                 q1 = 
            k1 /. a -> (k^2) + δ1 (e^2) + δ2 (e^4)  // 
              Collect[#, e, 
                  Simplify] &; 
                 q2 = {Coefficient[q1, e, 2] == 0,  
              Coefficient[q1, e, 4] == 
                0}; 
                 s1 = Solve[q2, {δ1, δ2}][[1]]; 
                 (k^2) + δ1 (e^2) + δ2 (e^4)  /. s1];

curves1=Table[ae1[j], {j, 0, 4}]

Extracting coefficients of   Cos[((k φ)/2)] in eq22 we obtain another infinite sequence of equations.

eq32 = (Coefficient[Coefficient[eq22 /. k -> 1, Sumk], 
              Cos[(φ/2)]] + 
            Coefficient[Coefficient[eq22 /. k -> 2, Sumk], 
              Cos[(φ/2)]]) == 0

coeffCos2[f_] := 
        Coefficient[f, 
            Cos[(φ/2) - k φ]] + (Coefficient[f, 
                Cos[((3  φ)/2) - k φ]]  /. 
              k -> k + 1) + (Coefficient[f, 
                Cos[(φ/2) + k φ]] /. 
              k -> k - 1) ;

eq33 = (coeffCos2[eq22]  /. Sumk -> 1  // 
            Collect[#, {c1[-3 + 2  k], c1[-1 + 2  k], c1[1 + 2  k]}, 
                Simplify] &) == 0

The corresponding matrix of coefficients can be written as

mat2[n_, eq0_Equal, eq1_Equal] := 
        Module[{j, l, 
            m}, 
                 m = {Table[Coefficient[eq0[[1]], 
                  c1[-1 + 2  j]], {j, 
                  n}]}; 
                 If[n > 1, Do[
              m = Append[m, 
                  Table[Coefficient[
                      eq1[[1]] /. 
                        k -> l, c1[-1 + 2  j]], {j, 
                      n}]], 
                              {l, 2, 
                n}]]; 
                           m ];

mat2[6, eq32, eq33]  // MatrixForm

The next recurrence relation is true for the determinants of the matrix above.

d2[1] = a - (1/4) + ((3  e)/8);
d2[
          2] = (a - (1/4) + ((3  e)/8)) (a - 
                (9/4)) + ((15  (e^2))/64);
d2[n_] := (a - ((((2  n - 1))^2)/4)) d2[n - 1] - 
          ((e^2)/64) (2  n + 1) (2  n - 1) (2  n - 3) (2  n - 
                5) d2[n - 2];

Using the recurrence relation we can easily find the stability boundary a = a(e) crossing the axis e = 0 in the point a = ((((2  k - 1))^2)/4). Again to find the corresponding function a = a(e) with an accuracy (e^(2  p)) we should use determinant of the (k+p)-th order. As a result we obtain the curves

ae2[k_] := 
        Module[{q1, q2, δ1, δ2, s1, k1},  
          k1 = d2[k + 
                2]; 
                 q1 = 
            k1 /. a -> ((((2  k - 1))^2)/4) + δ1 e + δ2 (e^2) + δ3 (e^3) + δ4 (e^4)  // 
              Collect[#, e, 
                  Simplify] &; 
                 q2 = {Coefficient[q1, e, 1] == 0, 
              Coefficient[q1, e, 2] == 0, 
              Coefficient[q1, e, 3] == 0,  
              Coefficient[q1, e, 4] == 
                0}; 
                 s1 = Solve[q2, {δ1, δ2, δ3, δ4}][[1]]; 
                 ((((2  k - 1))^2)/4) + δ1 e + δ2 (e^2) + δ3 (e^3) + δ4 (e^4)  /. s1];

curves2 = Table[ae2[j], {j, 4}]

The Odd Solution.

Substituting the fourth Fourier series in eq20 and separating the terms with even and odd k we have found

eq34 = 
        eq20 /. Z -> (Sumk (c1[2  k - 1] Sin[((2  k - 1)/2) #] + 
                        c1[2  k] Sin[k #]) &)  // 
          Simplify;
eq35 = reducing[eq34]

Extracting coefficients of  Sin[k φ]  in eq35 we obtain an infinite sequence of equations.

eq36 = (Coefficient[Coefficient[eq35 /. k -> 1, Sumk], 
              Sin[φ]] + 
            Coefficient[Coefficient[eq35 /. k -> 2, Sumk], 
              Sin[φ]]) == 0

coeffSin1[f_] := 
        Coefficient[f, 
              Sin[k φ]] - (Coefficient[f, 
                  Sin[φ - k φ]]  /. 
                k -> k + 1) + (Coefficient[f, 
                  Sin[φ + k φ]] /. 
                k -> k - 1) == 0 ;

eq37 = (coeffSin1[eq35]  /. Sumk -> 1  // 
          Collect[#, {c1[-2 + 2  k], c1[2  k], c1[2 + 2  k]}, 
              Simplify] &)

Then we define a function mat3 determining the matrix of the system.

mat3[n_, eq0_Equal, eq1_Equal] := 
        Module[{j, l, 
            m}, 
                 m = {Table[Coefficient[eq0[[1]], 
                  c1[2  j]], {j, 
                  n}]}; 
                 If[n > 1, Do[
              m = Append[m, 
                  Table[Coefficient[
                      eq1[[1]] /. 
                        k -> l, c1[2  j]], {j, 
                      n}]], 
                              {l, 2, 
                n}]]; 
                           m ];

The sixth order matrix for example can be written as

mat3[6, eq36, eq37]  // MatrixForm

Determinant of this matrix can be calculated with the next recurrence relation.

d3[1] = a - 1;
d3[2] = (a - 1) (a - 4);
d3[n_] := (a - (n^2)) d3[n - 1] - 
          ((e^2)/4) (n - 2) (n - 1) n (n + 1) d3[
              n - 2];

ae3[k_] := 
        Module[{q1, q2, δ1, δ2, s1, k1},  
          k1 = d3[k + 
                2]; 
                 q1 = 
            k1 /. a -> (k^2) + δ1 (e^2) + δ2 (e^4)  // 
              Collect[#, e, 
                  Simplify] &; 
                 q2 = {Coefficient[q1, e, 2] == 0,  
              Coefficient[q1, e, 4] == 
                0}; 
                 s1 = Solve[q2, {δ1, δ2}][[1]]; 
                 (k^2) + δ1 (e^2) + δ2 (e^4)  /. s1];

As a result w have obtained the curves a = a(e) in the form

curves3 = Table[ae3[j], {j, 4}]

Extracting coefficients of   Sin[((k φ)/2)]  in eq35 we again obtain an infinite sequence of equations.

eq38 = (Coefficient[Coefficient[eq35 /. k -> 1, Sumk], 
              Sin[(φ/2)]] + 
            Coefficient[Coefficient[eq35 /. k -> 2, Sumk], 
              Sin[(φ/2)]]) == 0

coeffSin2[
          f_] := -Coefficient[f, 
              Sin[(φ/2) - k φ]] - (Coefficient[f, 
                Sin[((3  φ)/2) - k φ]]  /. 
              k -> k + 1) + (Coefficient[f, 
                Sin[(φ/2) + k φ]] /. 
              k -> k - 1) ;

eq39 = (coeffSin2[eq35]  /. Sumk -> 1  // 
            Collect[#, {c1[-3 + 2  k], c1[-1 + 2  k], c1[1 + 2  k]}, 
                Simplify] &) == 0

Determining the matrix of the system and calculating its determinant we have obtained the last sequence of the curves a = a(e).

mat4[n_, eq0_Equal, eq1_Equal] := 
        Module[{j, l, 
            m}, 
                 m = {Table[Coefficient[eq0[[1]], 
                  c1[-1 + 2  j]], {j, 
                  n}]}; 
                 If[n > 1, Do[
              m = Append[m, 
                  Table[Coefficient[
                      eq1[[1]] /. 
                        k -> l, c1[-1 + 2  j]], {j, 
                      n}]], 
                              {l, 2, 
                n}]]; 
                           m ];

mat4[6, eq38, eq39]  // MatrixForm

Determinant of this matrix can be easily calculated with the next recurrence relation.

d4[1] = a - (1/4) - ((3  e)/8);
d4[
          2] = (a - (1/4) - ((3  e)/8)) (a - 
                (9/4)) + ((15  (e^2))/64);
d4[n_] := (a - ((((2  n - 1))^2)/4)) d4[n - 1] - 
          ((e^2)/64) (2  n + 1) (2  n - 1) (2  n - 3) (2  n - 
                5) d4[n - 2];

ae4[k_] := 
        Module[{q1, q2, δ1, δ2, s1, k1},  
          k1 = d4[k + 
                2]; 
                 q1 = 
            k1 /. a -> ((((2  k - 1))^2)/4) + δ1 e + δ2 (e^2) + δ3 (e^3) + δ4 (e^4)  // 
              Collect[#, e, 
                  Simplify] &; 
                 q2 = {Coefficient[q1, e, 1] == 0, 
              Coefficient[q1, e, 2] == 0, 
              Coefficient[q1, e, 3] == 0,  
              Coefficient[q1, e, 4] == 
                0}; 
                 s1 = Solve[q2, {δ1, δ2, δ3, δ4}][[1]]; 
                 ((((2  k - 1))^2)/4) + δ1 e + δ2 (e^2) + δ3 (e^3) + δ4 (e^4)  /. s1];

The corresponding sequence of the curves a = a(e) can be written as

curves4 = Table[ae4[j], {j, 4}]

Domains of Instability

As a result we have found the next stability boundaries.

curve = Union[curves1, curves2, curves3, curves4]

The corresponding curves are shown on the graph below.

p1 = 
        Plot[curve  // Evaluate, {e, 0, 1}, 
          PlotStyle -> Table[{th, Hue[0]}, {Length[curve]}], 
          PlotRange -> {{0, 0.9}, {0, 4.2}}, 
          AxesLabel -> {"e", "a"}];

We see that the domains of instability are only between the curves crossing the axis e = 0 in the points a = ((((2  k - 1))^2)/4). It is seen better from the graph below.

p2 = 
        FilledPlot[{curves2[[1]], curves4[[1]], curves2[[2]], 
            curves4[[2]]}, {e, 0, 1}, 
          Fills -> {{{1, 2}, Hue[0.5]}, {{3, 4}, Hue[0.5]}}, 
          Curves -> Back, AxesLabel -> {"e", "a"}, 
          DisplayFunction -> Identity];

p3 = 
        Plot[curve  // Evaluate, {e, 0, 1}, 
          PlotStyle -> Table[{th, Hue[0.5]}, {k, Length[curve]}], 
          DisplayFunction -> Identity];

Show[p2, p3, PlotRange -> {0, 4.5}, 
        DisplayFunction -> $DisplayFunction];

Stability of the solutions for perpendicular perturbations

Using the results above we can answer the question on stability of the solutions (4) with respect to the perpendicular perturbations. In order to do this let us define parameter a as a function of the number of the bodies n and the ratio μ = ((Subscript[m, 0])/m) taking into account that for every value of n there are n values of a. The corresponding function is denoted below as a1.

s1[
          n_] := (∑_(j = 1)^(n - 1))(1/Sin)[(π/n) j];
\(\(λ[n_, k_] := 
        2\ \(∑\+\(j = 1\)\%\(n - 1\)Sin[π\/n\ k\ j]\^2\/Sin[π\/n\ j]\^3\);\)\)
a1[n_, 
          k_, μ_] := ((λ[n, k] + 8  μ + 8 n If[k == n, 1, 0])/(2 (4  μ + s1[n])));

In the case of n = 6 for example we obtain the next list of parameters a as functions of μ.

p4[n_] := 
        Table[a1[n, k, μ], {k, n}];
p4[6]  // Simplify

The corresponding graphs are shown below.

Plot[p4[6] // Evaluate, {μ, 0, 20}, PlotRange -> All, 
        PlotStyle -> Table[{th, Hue[0.2  k]}, {k, Length[p4[6]]}], 
        AxesLabel -> {"μ", "a"}];

From the graphs above it's evident that parameters a = a(μ) are decreasing functions of μ which tends to the limit a = 1 for μ -> ∞. It means that for any value of the parameter e < 1 there is such value of (Subscript[μ, 0]) that for μ > (Subscript[μ, 0]) all values of a will be in the domain of stability on the ae-plane. Actually, let us first show the values of a = a(0)on the ae-plane.

p5 = 
        Plot[Table[p4[6][[
                  k]] /. μ -> 0, {k, 
                Length[p4[6]]}] // Evaluate, {e, 0, 1}, 
          PlotStyle -> Table[{th, Hue[0.2  k]}, {k, Length[p4[6]]}], 
          DisplayFunction -> Identity];

Show[p5, p2, p3, PlotRange -> {0, 4.5}, 
        AxesLabel -> {"e", "a"}, 
        DisplayFunction -> $DisplayFunction];

We see that there is an intersection of the line a = const and  a domain of instability. So in the case of μ = 0 for some values of parameter e a zero solution of the system (6) is unstable. On the next graph we show a domain of the function a = a(μ) changing in the case of μ ≥ 4.

p6 = 
        FilledPlot[{1, p4[
                  6][[6]] /. μ -> 4}, {e, 0, 1}, Fills -> {{{1, 2}, Hue[0.15]}}, 
          Curves -> Back, 
          DisplayFunction -> Identity];
Show[p6, p2, p3, PlotRange -> {0, 4.5}, 
        DisplayFunction -> $DisplayFunction];

Now we see that only for large enough values of the parameter e a zero solution of the system (6) is unstable. And with the parameter μ growth the domain of its stability will increase.

Conclusion

From the graphs above we can see that for large enough  μ a zero solution of equation eq12 is stable in linear approximation for any perturbations that are perpendicular to the plane of the configuration (Subscript[P, 0]) xy. It means that the central body of mass (Subscript[m, 0]) stabilizes the bodies motion with respect to the perpendicular perturbations.

IV. Studying Stability of the Solutions for the Perturbations Contained in the Plane of Orbits

To investigate linear stability of the solutions for the perturbations contained in the plane of particles' orbits it is convenient to rewrite equations eq10, eq11 in the normal form and introduce 4-dimensional vector

(Subscript[X, j])=(({{({{(Subscript[u, j])}, {(Subscript[γ, j])}})}, {({{(Subscript[ν, j])}, {(Subscript[ω, j])}})}}))

where

\(TextForm\`\(\(\ \)\(v\_j = \(d\ u\_j\)\/\(d\ φ\), \ \ ω\_j = \(\(\(d\ γ\_j\)\/\(d\ φ\)\)\(.\)\)\)\)\)

Then equations eq10, eq11 can be rewritten in matrix form as

((d (Subscript[X, j]))/(d φ)) = 
      A (φ) (Subscript[X, j]),

where the matrix A(φ) is a periodic  4 × 4  matrix of period 2π.

A (φ)=(({{({{0            , 0}, {0            , 0}}), ({{1, 0}, {0, 1}})}, {({{((Subscript[a, j])/(1 + e cosφ)), ((Subscript[b, j])/(1 + 
                            e cosφ))}, {-((Subscript[c, j])/(1 + 
                              e cosφ)), -((Subscript[d, j])/(1 + 
                              e cosφ))}}), ({{0, 2}, {-2, 0}})}}))

The problem of studying stability of a zero solution of the system (7) is very complicated and there is no any constructive method of its solving. Only in the case of e = 0 we have a linear system of differential equations with constant coefficients that is already investigated [4]. So the stability problem of a zero solution of the system (7) with periodic matrix is still open.

Acknowledgments

The author would like to thank Prof. Evgeni A. Grebenikov for very interesting and useful discussion of the stability problem.

References

1. B. Elmabsout. Sur l’existence de certaines configurations d’equilibre relatif dans le probleme des N corps. Celestial Mechanics, Vol. 41, (1988), 131-151. 
2. E.A.Grebenikov. The existence of the exact symmetric solutions in the plain Newton's problem of many bodies. Mathematical modeling, 10 (1998), No 8, 74-80.
3. E.A.Grebenikov. New exact solutions in the planar symmetrical (n+1)-body problem. Romanian Astronomical Journal, Vol. 7 (1997), No. 2, 151-156.  
4. B. Elmabsout. Stability of some degenerate positions of relative equilibrium in the n-body problem. Dynamics and Stability of Systems, Vol. 9 (1994), No. 4, 305-319.  
5. R.Grimshaw, Nonlinear Ordinary Diffrential Equations, CRC Press, 2000, 328 pp.

Generating Figures

Figure 1

n = 6;
e0 = 0;
α = ((2  π)/(n - 1));
r0[φ_] := 
        (1/(1 + e0 Cos[φ]));
x[i_] := 
        r0[π] Cos[α i + π];
y[i_] := 
        r0[π] Sin[α i + π];

\(\(obj = 
        Graphics[
          Join[{Line[
                Table[{x[j], y[j]}, {j, 0, 
                    n - 1}]]}, 
Table[{Hue[0.6], 
                PointSize[0.035], Point[{x[i], y[i]}]}, {i, 
                n - 1}], 
{Text[\*""(Subscript[P, 0])"", {0, \(-0.2\)}]}, {Text[\*""(Subscript[P, 1])"", {1.15  x[1], 
                  1.15  y[
                      1]}]}, 
{Text[\*""(Subscript[P, 2])"", {1.15  x[2], 1.15  y[2]}]}, {Text[\*""(Subscript[P, 3])"", {1.15  x[3], 
                  1.15  y[3]}]}, {Text[\*""(Subscript[P, 4])"", {1.15  x[4], 
                  1.15  y[4]}]}, {Text[\*""(Subscript[P, 5])"", {1.15  x[5], 
                  1.15  y[
                      5]}]}, 
{Text["Fig. 1", {0, \(-1.2\) r0[π]}]}, {PointSize[0.05], Hue[0], 
              Point[{0, 0}]}\ \ ]];\)\)

p1 = 
        Show[obj, AspectRatio -> Automatic, 
          DisplayFunction -> Identity];

p2 = 
        PolarPlot[
          Table[r0[φ + ((2  π)/(n - 1)) i], {i, n - 1}] // 
            Evaluate, {φ, 0, 2  π}, 
          PlotStyle -> Dashing[{0.03}], Axes -> None, 
          DisplayFunction -> Identity];

Show[p2, p1, DisplayFunction -> $DisplayFunction];

Figure 2

n = 6;
e0 = 0.7;
α = ((2  π)/(n - 1));
r0[φ_] := 
        (1/(1 + e0 Cos[φ]));
x[i_] := 
        r0[π] Cos[α i + π];
y[i_] := 
        r0[π] Sin[α i + π];

\(\(obj = 
        Graphics[
          Join[{Line[
                Table[{x[j], y[j]}, {j, 0, 
                    n - 1}]]}, 
Table[{Hue[0.6], 
                PointSize[0.035], Point[{x[i], y[i]}]}, {i, 
                n - 1}], 
{Text[\*""(Subscript[P, 0])"", {0, \(-0.4\)}]}, {Text[\*""(Subscript[P, 1])"", {1.15  x[1], 
                  1.15  y[
                      1]}]}, 
{Text[\*""(Subscript[P, 2])"", {1.15  x[2], 1.15  y[2]}]}, {Text[\*""(Subscript[P, 3])"", {1.15  x[3], 
                  1.15  y[3]}]}, {Text[\*""(Subscript[P, 4])"", {1.15  x[4], 
                  1.15  y[4]}]}, {Text[\*""(Subscript[P, 5])"", {1.15  x[5], 
                  1.15  y[
                      5]}]}, 
{Text["Fig. 2", {0, \(-1.2\) r0[π]}]}, {PointSize[0.05], Hue[0], 
              Point[{0, 0}]}\ \ ]];\)\)

p1 = 
        Show[obj, AspectRatio -> Automatic, 
          DisplayFunction -> Identity];

p2 = 
        PolarPlot[
          Table[r0[φ + ((2  π)/(n - 1)) i], {i, n - 1}] // 
            Evaluate, {φ, 0, 2  π}, 
          PlotStyle -> Dashing[{0.03}], Axes -> None, 
          DisplayFunction -> Identity];

Show[p2, p1, DisplayFunction -> $DisplayFunction];