Mathematica Day at the University of Warsaw
December 09, 2017

Application of the System Mathematica in the Sudy  
of the Circular Restricted Four-Body Problem

Alexander N. Prokopenya 
Warsaw University of Life Sciences - SGGW, Poland

E-mail: alexander_prokopenya@sggw.pl

Content

The problem statement

Equilibrium solutions

Stability analysis of the relative equilibrium positions

Normalization of the quadratic term (Subscript[H, 2])

Normalization of (Subscript[H, 3]) and third-order resonances

Normalization of (Subscript[H, 4]) and fourth-order resonances

The case f=(Subscript[c, 20])((Subscript[σ, 2])^2)+(Subscript[c, 11])(Subscript[σ, 1])(Subscript[σ, 2])+(Subscript[c, 02])((Subscript[σ, 1])^2)=0

Conclusion

Physical Problem

Three-body problem is a well-known model in Celestial Mechanics. It is known also that it is not integrable but there are five exact solutions of the problem which are known as the Lagrange triangle solutions and Euler collinear solutions. Three point particles (Subscript[P, 0]), (Subscript[P, 1]), (Subscript[P, 2]) having masses (Subscript[m, 0]), (Subscript[m, 1]), (Subscript[m, 2]), respectively, move uniformly on circular Keplerian orbits about their common center of mass. In the Lagrange case the particles are situated at the vertices of the equilateral triangle at any instant of time for any values of their masses. They attract each other according to the Newtonian law of gravitation and angular velocity of the particles is given by

((dν)/dt)≡(Subscript[ω, 0])=(((((G (Subscript[m, 0])(1+(Subscript[μ, 1])+(Subscript[μ, 2])))/((Subscript[r, 0])^3))))^(1/2)),  (Subscript[μ, 1])=((Subscript[m, 1])/(Subscript[m, 0])), (Subscript[μ, 2])=((Subscript[m, 2])/(Subscript[m, 0])),

where G is the constant of gravitation, ν is a polar angle, (Subscript[r, 0]) is a distance between the points (Subscript[P, 0]) and (Subscript[P, 1]).

In the Euler case all three particles also move uniformly on circular Keplerian orbits about their common center of mass but they are situated at the same straight line which rotates with constant angular velocity

((dν)/dt)≡(Subscript[ω, 0])=(((((G (Subscript[m, 0]))/((Subscript[r, 0])^3))(1+(Subscript[μ, 1])+(Subscript[μ, 2])((a/(|a((|)^3)))+((1-a)/(|1-a((|)^3)))))))^(1/2))

where parameter a is given by |a|=((|(Subscript[P, 0])(Subscript[P, 2])|)/(|(Subscript[P, 0])(Subscript[P, 1])|)).
For any values of these parameters there are three equilibrium positions of (Subscript[P, 2]) or three equilibrium values of the coordinate a:  a<0, 0<a<1, a>1 that are real roots of equation eq4.

Restricted problem of four bodies is formulated similarly to the restricted three-body problem. The fourth particle (Subscript[P, 3]) of negligible mass moves in the gravitational field of the primaries (Subscript[P, 0]), (Subscript[P, 1]), (Subscript[P, 2]), and the problem is to describe its motion. It is convenient to consider the particle (Subscript[P, 3]) motion in the rotating frame, where the particles (Subscript[P, 0]), (Subscript[P, 1]), (Subscript[P, 2]) are rest at the points (0,0), (1,0), (1/2, (Sqrt[3])/2) in the Lagrange case and at the pints  (0,0), (1,0), (a,0) in the Euler case, where there are three different cases: a<0, 0<a<1, a>1 and a=a((Subscript[μ, 1]),(Subscript[μ, 2])).

In case of the Lagrange triangle the massive particles  (Subscript[P, 0]), (Subscript[P, 1]), (Subscript[P, 2]) are situated at the vertices of equilateral triangle and their positions are the same for different values of parameters (Subscript[μ, 1]) and (Subscript[μ, 2]). 
Remind that there are eight, nine or ten equilibrium positions of the point (Subscript[P, 3]) in the plane Oxy, depending of parameters (Subscript[μ, 1]), (Subscript[μ, 2]). These positions are shown in picture below for (Subscript[μ, 1])=0.5 and 0≤(Subscript[μ, 2])≤(Subscript[μ, 1]).

Stability analysis has shown that equilibrium positions (Subscript[S, 2]), (Subscript[S, 3]), (Subscript[S, 5]),(Subscript[S, 6]), (Subscript[S, 8]) are unstable for any values of parameters (Subscript[μ, 1]), (Subscript[μ, 2]). The equilibrium positions (Subscript[S, 1]), (Subscript[S, 4])  and (Subscript[S, 7]) are stable in Liapunov's sense for all values of  the parameters (Subscript[μ, 1]), (Subscript[μ, 2]) from the domain of its linear stability except for 
a) the curve (Subscript[σ, 1])=2(Subscript[σ, 2]), where the third order resonance takes place;
b) the curve (Subscript[σ, 1])=3(Subscript[σ, 2]), where the fourth order resonance takes place (equilibrium position (Subscript[S, 4]) is stable in this case);
c) two points in the plane (Subscript[μ, 1])(Subscript[Oμ, 2]), where the curve f=0 crosses the resonance curves (Subscript[σ, 1])=4(Subscript[σ, 2]) and (Subscript[σ, 1])=5(Subscript[σ, 2]), respectively.

In the Euler case the massive particles  (Subscript[P, 0]), (Subscript[P, 1]), (Subscript[P, 2]) are situated at the axis Ox and position of the particle (Subscript[P, 2]) depends on the values of parameters (Subscript[μ, 1]) and (Subscript[μ, 2]). For any values of these parameters there are three equilibrium positions of (Subscript[P, 2]) or three equilibrium values of the coordinate a:  a<0, 0<a<1, a>1 that are real roots of equation eq4.

eq4=(((-1+a))^3) (a^4) (1+(Subscript[μ, 1])+(((1-a)/((((((-1+a))^2)))^(3/2)))+(a/((((a^2)))^(3/2)))) (Subscript[μ, 2]))==(a^3) (-1+(Sqrt[(((-1+a))^2)])+3 a-3 (a^2)+(a^3)) (Subscript[μ, 1])+(((-1+a))^3) (Sqrt[(a^2)]) (1+(Subscript[μ, 2]));

In the domain a>1 equilibrium position of (Subscript[P, 2])  is determined as a root of equation eq4a.

eq4a=eq4 //Simplify[#,a>1]&;
sol4a[μ1_,μ2_]:=a/.FindRoot[eq4a /.{(Subscript[μ, 1])->μ1, (Subscript[μ, 2])->μ2},{a,2}]

Manipulate[Plot[eq4a[[1]]-eq4a[[2]] /.{(Subscript[μ, 1])->μ1,(Subscript[μ, 2])->μ2},{a,1,2},PlotRange->6,BaseStyle->{FontSize->14}],{μ1,1/100,10},{μ2,1/1000,1}]

In the domain 0<a<1 equilibrium position of (Subscript[P, 2])  is determined as a root of equation eq4b.

eq4b=eq4 //Simplify[#,0<a<1]&;
sol4b[μ1_,μ2_]:=a/.FindRoot[eq4b /.{(Subscript[μ, 1])->μ1, (Subscript[μ, 2])->μ2},{a,1/2}]

Manipulate[Plot[eq4b[[1]]-eq4b[[2]] /.{(Subscript[μ, 1])->μ1,(Subscript[μ, 2])->μ2},{a,0,1},PlotRange->1,BaseStyle->{FontSize->14}],{μ1,1/100,1},{μ2,1/100,1}]

In the domain a<0 equilibrium position of (Subscript[P, 2])  is determined as a root of equation eq4c.

eq4c=eq4 //Simplify[#,a<0]&;
sol4c[μ1_,μ2_]:=a/.FindRoot[eq4c /.{(Subscript[μ, 1])->μ1, (Subscript[μ, 2])->μ2},{a,-1}]

Manipulate[Plot[eq4c[[1]]-eq4c[[2]] /.{(Subscript[μ, 1])->μ1,(Subscript[μ, 2])->μ2},{a,-1,0},PlotRange->5,BaseStyle->{FontSize->14}],{μ1,1/100,10},{μ2,1/100,1}]

Using dimensionless coordinates and momenta and considering the polar angle ν as an independent variable, one can write the Hamiltonian of the system in the form

H=(((Px[ν])^2)/2)+(((Py[ν])^2)/2)-Py[ν] x[ν]+Px[ν] y[ν]-(1/(κ (Sqrt[((x[ν])^2)+((y[ν])^2)])))+((Subscript[μ, 1])/κ) (x[ν]-(1/(Sqrt[1-2 x[ν]+((x[ν])^2)+((y[ν])^2)])))+((Subscript[μ, 2])/κ) ((((Sqrt[(a^2)]) x[ν])/(a^3))-(1/(Sqrt[(((a-x[ν]))^2)+((y[ν])^2)])));
solκ=κ->1+(Subscript[μ, 1])+(( (Subscript[μ, 2]))/(a (1-a)))((a/(Sqrt[(((-1+a))^2)]))+(1/(Sqrt[(a^2)]))-((Sqrt[(a^2)])/a));

Then equations of motion are obtained in a standard way.

eq6={x'[ν]==D[H,Px[ν]],y'[ν]==D[H,Py[ν]], 
    Px'[ν]==-D[H,x[ν]],Py'[ν]==-D[H,y[ν]]} //Simplify;

Equilibrium solutions

We are looking for stationary solutions of the equations of motion.

eq7=eq6/.{x->(x&),y->(y&),Px->(px&),Py->(py&)};
eq8=Simplify[eq7[[{3,4}]]/.Solve[eq7[[{1,2}]],{px,py}][[1]]];
eq9=Thread[Simplify[Expand[(Thread[κ (#1[[1]]-#1[[2]])]&)/@eq8]]==0]

Equations eq9 determine two curves in the Oxy plane and the points of their intersection corresponds to the equilibrium solutions. To estimate the number of equilibrium points and their positions in the Oxy plane we plot two curves using ContourPlot function.

For small values of parameter (Subscript[μ, 2]) each point of libration (Subscript[L, 1]), (Subscript[L, 2]), (Subscript[L, 3]) generates two equilibrium solutions of the particle  (Subscript[P, 3]) if the particle (Subscript[P, 2]) is situated in the corresponding libration point (in the figure above the points (Subscript[S, 2]) and (Subscript[S, 3]) arise from the point (Subscript[L, 2])).

Manipulate[ContourPlot[eq9 /.solκ/.{(Subscript[μ, 1])->μ1,(Subscript[μ, 2])->μ2} /. a->sol4c[μ1,μ2] //Evaluate,{x,-2,2.},{y,-2,2},ContourStyle->{{Thick,Red},{Thick,Dashed,Blue}},Epilog->{{PointSize[0.04],Point[{0,0}]},{PointSize[0.03],Point[{1,0}]},{PointSize[0.02],Point[{sol4c[μ1,μ2] ,0}]},
Text[StringForm["(Subscript[μ, 1]) = ``",μ1],{-1,1.55}],Text[StringForm["(Subscript[μ, 2]) = ``",μ2],{-1,1.35}]},BaseStyle->{FontSize->14}],{μ1,0.00001,1},{μ2,0.00001,1}]

In general, there are six equilibrium positions for any values of the parameters (Subscript[μ, 1]), (Subscript[μ, 2]).

In the case of (Subscript[μ, 2])=0 there exist five equilibrium positions what is known from the restricted three-body problem. For small values of (Subscript[μ, 2]) and some fixed value of (Subscript[μ, 1]) the equilibrium positions can be found in the form of power series in (Subscript[μ, 2]). The following procedure finds the corresponding power series accurate up to n-th order, where (x0,y0) is the starting point.

Function curves16[μ1,n,x0,y0,sola] for constructing a single solution starting at the point ((Subscript[x, 0]),(Subscript[y, 0])) for different values of (Subscript[μ, 1]) and different values of a at (Subscript[μ, 2])=0 (sola determines one of the three possible initial positions of the body (Subscript[P, 2]): sol4a[(Subscript[μ, 1]),0] corresponds to the libration point (Subscript[L, 2]), sol4b[(Subscript[μ, 1]),0] corresponds to the libration point (Subscript[L, 1]), sol4c[(Subscript[μ, 1]),0] corresponds to the libration point (Subscript[L, 3])) in the form of a polynomial in term of (Subscript[μ, 2]) of degree n.

Three solutions starting from the libration point (Subscript[L, 4]) when (Subscript[P, 2]) is situated in one of the points (Subscript[L, 1]), (Subscript[L, 2]), (Subscript[L, 3]) are

{curves16[1/5,3, 1/2,(Sqrt[3])/2,sol4a[1/5,0]],curves16[1/5,3, 1/2,(Sqrt[3])/2,sol4b[1/5,0]],curves16[1/5,3, 1/2,(Sqrt[3])/2,sol4c[1/5,0]]}//Column

Two collinear solutions starting from the libration points (Subscript[L, 1]) and (Subscript[L, 3]) when (Subscript[P, 2]) is situated at the point  (Subscript[L, 2])

{curves16[1/5,3, sol4b[1/5,0],0,sol4a[1/5,0]], curves16[1/5,3, sol4c[1/5,0],0,sol4a[1/5,0]]}//Column

Function curvesMul[μ1,n,sola] for constructing multiple solutions starting from the point (a,0) for (Subscript[μ, 2])=0, where the body (Subscript[P, 2]) is situated. There are three such points at the axis Ox, they are the points of libration (Subscript[L, 1]), (Subscript[L, 2]), (Subscript[L, 3]) for (Subscript[μ, 2])=0.

One can readily see that two equilibrium points arise from the points of libration (Subscript[L, 1]), (Subscript[L, 2]), (Subscript[L, 3]) for (Subscript[μ, 2])>0.

curvesMul[1/5,2,sol4c[1/5,0]]

Positions of the particles (Subscript[P, 2]) and (Subscript[P, 3]) for (Subscript[μ, 1])=0.2 and different values of 0≤(Subscript[μ, 2])≤1

Figure below shows positions of the body (Subscript[P, 3]) for different values of parameter (Subscript[μ, 2]) and (Subscript[μ, 1])=0.2 when the body (Subscript[P, 2]) is situated in one of the libration points (Subscript[L, 1]), (Subscript[L, 2]), (Subscript[L, 3]).

Stability analysis of equilibrium positions

Arnold-Moser Theorem

H

Let the Hamiltonian be represented in the form

H=(Subscript[σ, 1])(Subscript[r, 1])-(Subscript[σ, 2])(Subscript[r, 2])+(Subscript[c, 20])((Subscript[r, 1])^2)+(Subscript[c, 11])(Subscript[r, 1])(Subscript[r, 2])+(Subscript[c, 02])((Subscript[r, 2])^2)+(Subscript[H, 5])((Subscript[r, j]),(Subscript[φ, j]))+...,

where (Subscript[σ, 1]),(Subscript[σ, 2]) are real numbers which are frequencies for the linearized dynamical system, and the following conditions are satisfied:

a)   (Subscript[k, 1]) (Subscript[σ, 1])+(Subscript[k, 2]) (Subscript[σ, 2])≠0  for all integers (Subscript[k, 1]), (Subscript[k, 2]), |(Subscript[k, 1])|+|(Subscript[k, 2])|≤4;

b)   f=(Subscript[c, 20])((Subscript[σ, 2])^2)+(Subscript[c, 11]) (Subscript[σ, 1]) (Subscript[σ, 2])+(Subscript[c, 02]) ((Subscript[σ, 1])^2)≠0.

Then equilibrium solution of the Hamiltonian system is stable in the Liapunov sense.

Third-order resonance (A.P. Markeev)

Let the Hamiltonian of the system

H=(Subscript[H, 2])+(Subscript[H, 3])+...,

be represented in the form

H=(Subscript[σ, 1])(Subscript[r, 1])-(Subscript[σ, 2])(Subscript[r, 2])+B Cos[(Subscript[φ, 1])+2 (Subscript[φ, 2])] (Sqrt[(Subscript[r, 1])]) (Subscript[r, 2])+(Subscript[H, 4])((Subscript[r, j]),(Subscript[φ, j]))+...,

where (Subscript[σ, 1])=2(Subscript[σ, 2]) or 2(Subscript[σ, 1])=(Subscript[σ, 2]).

Then equilibrium solution of the Hamiltonian system is unstable if B≠0.
If B=0 but (Subscript[c, 20])+2(Subscript[c, 11])+4(Subscript[c, 02])≠0 then equilibrium solution is stable in the Liapunov sense.

Fourth-order resonance (A.P. Markeev)

Let the Hamiltonian be represented in the form

H=(Subscript[σ, 1])(Subscript[r, 1])-(Subscript[σ, 2])(Subscript[r, 2])+(Subscript[c, 20])((Subscript[r, 1])^2)+(Subscript[c, 11])(Subscript[r, 1])(Subscript[r, 2])+(Subscript[c, 02])((Subscript[r, 2])^2)+A Cos[(Subscript[φ, 1])+3 (Subscript[φ, 2])] (Sqrt[(Subscript[r, 1])]) ((Subscript[r, 2])^(3/2))+O(((((Subscript[r, 1])+(Subscript[r, 2])))^(5/2))),

where (Subscript[σ, 1])=3(Subscript[σ, 2]).

Then equilibrium solution of the Hamiltonian system is unstable for |(Subscript[c, 20])+3(Subscript[c, 11])+9(Subscript[c, 02])|<3(Sqrt[3])|A|
and is stable in the Liapunov sense for |(Subscript[c, 20])+3(Subscript[c, 11])+9(Subscript[c, 02])|>3(Sqrt[3])|A|.

Normalization of the Hamiltonian

Let us add small perturbation to the equilibrium solution.

rul2={x[ν]->(x0+δ x[ν]),y[ν]->(y0+δ y[ν]), Px[ν]->(-y0+δ Px[ν]), Py[ν]->(x0+δ Py[ν])};

Then we can expand the Hamiltonian into power series in the neighbourhood of the equilibrium solution

Hexp=Series[H /.rul2,{δ,0,4}] //Normal;

One can readily see that zero order term (Subscript[H, 0]) is a constant which can be neglected because it doesn't influence the equations of motion.

H0=Coefficient[Hexp,δ,0] //Simplify

Coefficient of δ determines the first order term (Subscript[H, 1]) that is equal to zero due to equilibrium condition eq9.

H1=Coefficient[Hexp,δ,1] //Collect[#,{x[ν],y[ν]},Simplify]&

Therefore, the first nonzero term in the Hamiltonian expansion is the quadratic one (Subscript[H, 2]) that is given by

H2=Coefficient[Hexp,δ,2] //Collect[#,{Px[ν],Py[ν],x[ν],y[ν]},Simplify]&

To demonstrate main steps in normalizing the Hamiltonian and to simplify the calculations we'll use the following notations

H20=(((Px[ν])^2)/2)+(((Py[ν])^2)/2)-Py[ν] x[ν]+Px[ν] y[ν]+h20  ((x[ν])^2)+h11 x[ν] y[ν]+h02  ((y[ν])^2);

H30=h330  ((x[ν])^3)+h321 ((x[ν])^2) y[ν]+h312 x[ν] ((y[ν])^2)+h303  ((y[ν])^3);

H40=h440  ((x[ν])^4)+h431 ((x[ν])^3) y[ν]+h422 ((x[ν])^2) ((y[ν])^2)+h413 x[ν] ((y[ν])^3)+h404  ((y[ν])^4);

Characteristic exponents

To calculate characteristic exponents of the system let us rewrite the linearized equations in the vector form.

X2={x[ν], y[ν],Px[ν],Py[ν]};
J2={{0,0,1,0},{0,0,0,1},{-1,0,0,0},{0,-1,0,0}};
H2mat=Table[D[H20,X2[[j]],X2[[k]] ],{j,4},{k,4}];
D[X2,ν]==J2.H2mat.X2 //Thread

Obviously, the characteristic exponents are the eigenvalues of the matrix J.(Subscript[H, 2mat]) and can be found as zeros of the characteristic polynomial.

pol=CharacteristicPolynomial[J2.H2mat,λ] //Collect[#,λ]&

We are looking for purely imaginary characteristic exponents (Subscript[λ, 1,3])=±i (Subscript[σ, 1]), (Subscript[λ, 2,4])=±i (Subscript[σ, 2]), where two frequencies (Subscript[σ, 1]), (Subscript[σ, 2]) are given by

eigenval1=Solve[(pol /.λ->I σ)==0,σ] //Simplify;
eigenval={(Subscript[σ, 1])->(σ /.eigenval1[[4]]), (Subscript[σ, 2])->(σ /.eigenval1[[2]])}

The system has four different purely imaginary characteristic exponents if the following conditions are satisfied
1-2 h02-(h11^2)-2 h20+4 h20 h02>0 ,  
1+ h20+h02>0
4 h20+(h20^2)+(h02^2)+4 h02-2 h20 h02+(h11^2)>0

Analyzing these inequalities, we have found that at least one of them is not satisfied for the collinear points (Subscript[S, 1]), (Subscript[S, 2]), (Subscript[S, 3]), (Subscript[S, 4]) and, hence, these four equilibrium points are unstable.

In case of the point (Subscript[S, 5]) (for (Subscript[S, 6])) it is the same there is a domain of linear stability bounded by the curve (Subscript[σ, 1])=(Subscript[σ, 2]), it is shown at the figure below (dashed curve) together with four resonance curves when the particle (Subscript[P, 2]) is situated in the libration points (Subscript[L, 1]) or (Subscript[L, 2]).

Below the domain of stability bounded by the curve  (Subscript[σ, 1])=(Subscript[σ, 2]) (dashed curve) is shown in the case when the particle (Subscript[P, 2]) is situated in the libration points (Subscript[L, 3]). One can readily see that this figure is symmetrical with respect to the line (Subscript[μ, 2])=(Subscript[μ, 1]).

Normalization of the quadratic term (Subscript[H, 2])

Normalization of the quadratic term (Subscript[H, 2]) means that we are looking for the real-valued canonical transformation (Subscript[X, 2])=A.X, (X^T)=((Subscript[q, 1]),(Subscript[q, 2]),(Subscript[p, 1]),(Subscript[p, 2])), reducing the second order polynomial (Subscript[H, 2]) to the form

(Subscript[H, 2])=±((Subscript[σ, 1])/2)(((Subscript[p, 1])^2)+((Subscript[q, 1])^2))±((Subscript[σ, 2])/2)(((Subscript[p, 2])^2)+((Subscript[q, 2])^2))

where the sign "+" or "–" of each term is determined as a result of normalization.
Matrix A of such transformation can be constructed with a usage of the eigenvectors of the matrix J.(Subscript[H, 2mat]).

X={(Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4])};
eqEigen=J2.H2mat.X == λ X //Thread;
solEigen=Solve[eqEigen[[{1,2,3} ]],X[[{2,3,4}]]][[1]] //Simplify;
vectors=X  /.solEigen /.(Subscript[x, 1])->1

Obviously, the vector above determines two pairs of complex conjugate eigenvectors corresponding to the eigenvalues (Subscript[λ, 1])=-(Subscript[λ, 3])=i (Subscript[σ, 1]), (Subscript[λ, 2])=-(Subscript[λ, 4])=i (Subscript[σ, 2]).

To construct the matrix A of canonical transformation we substitute λ->i (Subscript[σ, k])  (k=1,2) into the general expression for the eigenvector vectors multiplied by a real constant (Subscript[c, k]) (k=1,2) and extract real and imaginary parts of each eigenvector.

re[k_]=c[k]ComplexExpand[Re[vectors /.λ->I (Subscript[σ, k])]] //Simplify

im[k_]=c[k]ComplexExpand[Im[vectors /.λ->I (Subscript[σ, k])]] //Simplify

Unknown coefficients c[k] are found from the normalization condition.

sol2=Solve[((re[k].J2.im[k] //Simplify) /.((c[k])^2)->cc)==(1/4),cc ][[1]];
sol2 /. cc->((c[k])^2)

As coefficients (Subscript[c, k]) are real numbers their square ((c[k])^2) must be positive. Calculation shows that ((c[1])^2)>0 but ((c[2])^2)<0  for any values of (Subscript[μ, 1]), (Subscript[μ, 2]) what means that c[2] is an imaginary number.

Making the replacement (Subscript[σ, 2])->-(Subscript[σ, 2]) (it means we should substitute (Subscript[λ, 2])=-i (Subscript[σ, 2]) in the expression for the eigenvector instead of (Subscript[λ, 2])=i (Subscript[σ, 2])), we obtain  ((c[2])^2)>0.

norm={c[1]->((Sqrt[cc]) /.sol2 /.k->1),c[2]->((Sqrt[cc]) /.sol2 /.k->2 /.(Subscript[σ, 2])->-(Subscript[σ, 2]))} //Simplify;

Then the matrix A is determined as

AA=Transpose[{-2im[1],-2im[2],2re[1],2re[2]}]/.(Subscript[σ, 2])->-(Subscript[σ, 2])   //Simplify

One can easily check that this matrix is symplectic and, hence, the transformation (Subscript[X, 2])=A.X, (X^T)=((Subscript[q, 1]),(Subscript[q, 2]),(Subscript[p, 1]),(Subscript[p, 2])) is canonical.

Transpose[AA].J2.AA  /.norm /.eigenval//Simplify

The transformed matrix (SuperStar[H])=-J (A^(-1))J H A is diagonal.

(SuperStar[H])=-J (A^(-1))J H A=(({{(Subscript[σ, 1]), 0, 0, 0}, {0, -(Subscript[σ, 2]), 0, 0}, {0, 0, (Subscript[σ, 1]), 0}, {0, 0, 0, -(Subscript[σ, 2])}}))

Actually, doing the calculations, we obtain

H2norm=-J2.Inverse[AA].J2.H2mat.AA /.eigenval//FullSimplify

Then the second order term in the Hamiltonian takes the normalized form

H2norm=(1/2) (Subscript[σ, 1]) (p1^2+q1^2)-(1/2) (Subscript[σ, 2]) (p2^2+q2^2);

Note that we have obtained the sign "–" for the second term as a result of quadratic normalization. 
One can readily see that (Subscript[H, 2]) is neither positive nor negative definite function. Therefore, we cannot conclude on stability or instability of equilibrium position using the principle of linearized stability.
So we have to solve the stability problem in a strict nonlinear formulation, using Arnold-Moser and Markeev theorems.

Normalization of (Subscript[H, 3])

On substituting the canonical transformation (Subscript[X, 2])=A.X normalizing the term (Subscript[H, 2]) into the third order term (Subscript[H, 3]) it takes the form

H31=Sum[If[i+j+k+l==3,h3[i,j,k,l] (q1^i) (q2^j) (p1^k) (p2^l),0],{i,0,3},{j,0,3},{k,0,3},{l,0,3}];

Normalization of (Subscript[H, 3]) means that we construct such canonical transformation that (Subscript[H, 3]) would be cancelled. Following to Birkhoff, we seek the corresponding generating function in the form of the third order polynomial

S3[p1n_,p2n_,q1_,q2_]=Sum[If[i+j+k+l==3,s3[i,j,k,l] (q1^i) (q2^j) (p1n^k) (p2n^l),0],{i,0,3},{j,0,3},{k,0,3},{l,0,3}];

where the coefficients s3[i,j,k,l] are unknown. Then old canonical variables q1, q2, p1, p2 and new ones q1n, q2n, p1n, p2n are connected by means of the relationships

eq3={q1n==q1+D[S3[p1n,p2n,q1,q2],p1n], 
      q2n==q2+D[S3[p1n,p2n,q1,q2],p2n], 
      p1==p1n+D[S3[p1n,p2n,q1,q2],q1], 
      p2==p2n+D[S3[p1n,p2n,q1,q2],q2] } ;

Solving eq3, we obtain old variables in the form of power series in new variables. To normalize (Subscript[H, 3]) it is sufficient to find the corresponding solution up to quadratic terms. The corresponding canonical transformation is given by

rul6={q1->q1-D[S3[p1,p2,q1,q2],p1],
           q2->q2-D[S3[p1,p2,q1,q2],p2],
           p1->p1+D[S3[p1,p2,q1,q2],q1],
           p2->p2+D[S3[p1,p2,q1,q2],q2]};

Now we can substitute rul6 into H2norm and H3norm and calculate the Taylor series expansion of the Hamiltonian in powers of q1, q2, p1, p2. One can readily check that the second order term H2new in this expansion will have the some form as H2norm form while the third order term in H23new will be again a sum of twenty terms of the form  h3[i,j,k,l] (q1^i) (q2^j) (p1^k) (p2^l).

rul7={p1->δ p1,p2->δ p2,q1->δ q1,q2->δ q2};
H23new=H2norm+H31 /.rul6/.rul7 //Series[#,{δ,0,3}]& //Normal ;

Now we can extract new coefficients h3[i,j,k,l] in the expansion H23new.

var5={};Do[If[i+j+k+l==3,var5=Append[var5,(q1^i) (q2^j) (p1^k) (p2^l)]],{i,0,3},{j,0,3},{k,0,3},{l,0,3}];
eq6=Table[Coefficient[H23new /.δ->1,var5[[j]] ],{j,Length[var5]}]; eq6 //Simplify//Column

These 20 expressions contain 20 unknown coefficients  s3[i,j,k,l] which determine the generating function (Subscript[S, 3])((Subscript[p, 1]),(Subscript[p, 2]),(Subscript[q, 1]),(Subscript[q, 2])). A list of these coefficients is shown below.

var6={};Do[If[i+j+k+l==3,var6=Append[var6,s3[i,j,k,l]]],{i,0,3},{j,0,3},{k,0,3},{l,0,3}];var6

Note that we have to find such values for the coefficients s3[i,j,k,l] which eliminate the third-order term in the (Subscript[H, 3]). It means that all coefficients h3[i,j,k,l] in H23new must be equal to zero.
Therefore, we can try to solve the system of twenty equations of the form  h3[i,j,k,l]=0 and to find the corresponding coefficients s3[i,j,k,l].

Analysis of the coefficients h3[i,j,k,l] shows that in fact we have five independent subsystems for determination of the coefficients s3[i,j,k,l]. The first sybsystem is formed by three coefficients of h3[0,0,1,2], h3[0,2,1,0], h3[1,1,0,1]. It determines the coefficients s3[0,1,1,1], s3[1,0,0,2], s3[1,2,0,0] and is given by

eq6a=eq6[[{2,9,14}]];eq6a //Simplify//Column

Note that variables s3[0,1,1,1], s3[1,0,0,2], s3[1,2,0,0] appear only in the expressions for coefficients h3[0,0,1,2], h3[0,2,1,0], h3[1,1,0,1]. Equating these coefficients to zero and solving the system obtained, we find

sol6a=Solve[(eq6a==0 //Thread),{s3[0,1,1,1],s3[1,0,0,2],s3[1,2,0,0]}][[1]] //Simplify

One can readily see that this solution exists only if (Subscript[σ, 1])≠0 and (Subscript[σ, 1])±2(Subscript[σ, 2])≠0.

Analyzing the rest four subsystems of eq6, we have shown that if the conditions

(Subscript[σ, 1])≠0, (Subscript[σ, 2])≠0,  (Subscript[σ, 1])±2(Subscript[σ, 2])≠0,  2(Subscript[σ, 1])±(Subscript[σ, 2])≠0

are satisfied, what means that there are no resonances of the first and third order in the system, then the third order term (Subscript[H, 3]) in the Hamiltonian expansion is eliminated by means of the canonical transformation determined by the generating function S3 we have constructed above.

Note that in the considered restricted four-body problem there are such values of parameters (Subscript[μ, 1]) and (Subscript[μ, 2]) for which the condition of third-order resonance (Subscript[σ, 1])-2(Subscript[σ, 2])=0  is satisfied. In such a case we cannot equate coefficients  h3[0,0,1,2], h3[0,2,1,0], h3[1,1,0,1] to zero and eliminate the corresponding terms in (Subscript[H, 3]). 
Analyzing the system eq6a one can show that its determinant is equal to zero but there are nonzero minors of the second order. Therefore, two variables of s3[0,1,1,1], s3[1,0,0,2], s3[1,2,0,0] are independent and we can introduce one more variable (Subscript[B, 1])according to the following rule

eq6a1=eq6a=={(B1/(2(Sqrt[2]))),-(B1/(2(Sqrt[2]))),-(B1/(Sqrt[2]))} //Thread; eq6a1

We make such a choice to be able to represent a final expression for the Hamiltonian in the form which accepts application of A.P. Markeev's theorem on stability of equilibrium position in case of the third order resonance.

Then we solve the system eq6a1 with respect to s3[0,1,1,1], s3[1,0,0,2], for example, and B1

sol6a1=Solve[eq6a1 /.(Subscript[σ, 1])->2(Subscript[σ, 2]),{s3[1,0,0,2],s3[0,1,1,1],B1}][[1]]//Simplify

and the third variable s3[1,2,0,0] remains as free parameter.

The same calculations are done for the second subsystem whose determinant equals to zero under the third order resonance (Subscript[σ, 1])-2(Subscript[σ, 2])=0.

eq6b=eq6[[{6,11,16}]]; eq6b1=eq6b=={(B2/(Sqrt[2])),(B2/(2(Sqrt[2]))),-(B2/(2(Sqrt[2])))} //Thread;
sol6b1=Solve[eq6b1 /.(Subscript[σ, 1])->2(Subscript[σ, 2]),{s3[0,0,1,2],s3[1,1,0,1],B2}][[1]]//Simplify

The rest 14 coefficients h3[i,j,k,l] become zeros after normalization of (Subscript[H, 3]) and after some standard calculations the Hamiltonian takes the form

H=2 (Subscript[r, 1]) (Subscript[σ, 2])-(Subscript[r, 2]) (Subscript[σ, 2])+(Sqrt[B1^2+B2^2]) (Sqrt[(Subscript[r, 1])]) (Subscript[r, 2]) Cos[(Subscript[φ, 1])+2 (Subscript[φ, 2])]+(Subscript[H, 4])+...

where (Subscript[r, 1]), (Subscript[r, 2]), (Subscript[φ, 1]), (Subscript[φ, 2]) are new canonical variables of the "action-angle" type. According to Markeev's theorem, an equilibrium position is unstable if (Sqrt[((Subscript[B, 1])^2)+((Subscript[B, 2])^2)])≠0.

We checked this condition for the point (Subscript[S, 5]) and three equilibrium positions (Subscript[L, 1]), (Subscript[L, 2]), (Subscript[L, 3]) of the particle (Subscript[P, 2]) that is stable in linear approximation and showed that in all three cases it is unstable under the third order resonance of the form (Subscript[σ, 1])-2(Subscript[σ, 2])=0.

Normalization of (Subscript[H, 4])

It is assumed that there are no resonances of the first and third order and normalization of the second and third order terms (Subscript[H, 2]), (Subscript[H, 3]) have been done. Then second-order term (Subscript[H, 2]) has the form

H2norm=(1/2) (Subscript[σ, 1]) (p1^2+q1^2)-(1/2) (Subscript[σ, 2]) (p2^2+q2^2);

third-order term (Subscript[H, 3])=0 and the fourth order term (Subscript[H, 4]) is represented in a general form

H41=Sum[If[i+j+k+l==4,h4[i,j,k,l] (q1^i) (q2^j) (p1^k) (p2^l),0],{i,0,4},{j,0,4},{k,0,4},{l,0,4}]

where 35 coefficients h4[i,j,k,l] are determined by the parameters (Subscript[μ, 1]), (Subscript[μ, 2]) and coefficients s3[i,j,k,l] of the generating function (Subscript[S, 3]) we have found above.

Normalization of (Subscript[H, 4]) is done similarly to the previous case of the normalization of (Subscript[H, 3]). The main difference is that now the generating function of the corresponding canonical transformation is constructed in the form of fourth-order polynomial

S4[p1_,p2_,q1_,q2_]=Sum[If[i+j+k+l==4,s4[i,j,k,l] (q1^i) (q2^j) (p1^k) (p2^l),0],{i,0,4},{j,0,4},{k,0,4},{l,0,4}];

where the coefficients s4[i,j,k,l] are unknown. Then old canonical variables q1, q2, p1, p2 and new ones q1n, q2n, p1n, p2n are connected by means of the relationships

eq7={q1n==q1+D[S4[p1n,p2n,q1,q2],p1n],  
       q2n==q2+D[S4[p1n,p2n,q1,q2],p2n], 
         p1==p1n+D[S4[p1n,p2n,q1,q2],q1],
         p2==p2n+D[S4[p1n,p2n,q1,q2],q2] } ;

Solving eq7, we obtain old variables in the form of power series in new variables. To normalize (Subscript[H, 4]) it is sufficient to find the corresponding solution up to cubic terms and to obtain new coefficients h4[i,j,k,l].

Analysis of the coefficients h4[i,j,k,l] shows that in the absence of resonances up to the fourth order inclusively the Hamiltonian is reduced to the form

H=(Subscript[r, 1]) (Subscript[σ, 1])-(Subscript[r, 2]) (Subscript[σ, 2])+(Subscript[c, 20]) (Subscript[r, 1]^2)+(Subscript[c, 02]) (Subscript[r, 2]^2)+(Subscript[c, 11]) (Subscript[r, 1]) (Subscript[r, 2])+(Subscript[H, 5])+...

where the coefficients (Subscript[c, 20]), (Subscript[c, 11]), (Subscript[c, 02]) are very cumbersome and so we do not show them here. Such a form of the Hamiltonian accepts applying the Arnold-Moser theorem.

In the case of fourth order resonance of the form (Subscript[σ, 1])-3(Subscript[σ, 2])=0 which takes place in the system under consideration the Hamiltonian is reduced to the form

H=3(Subscript[r, 1]) (Subscript[σ, 2])-(Subscript[r, 2]) (Subscript[σ, 2])+(Subscript[c, 20]) (Subscript[r, 1]^2)+(Subscript[c, 02]) (Subscript[r, 2]^2)+(Subscript[c, 11]) (Subscript[r, 1]) (Subscript[r, 2])+(Subscript[r, 2])(Sqrt[(Subscript[r, 1])(Subscript[r, 2])])(Sqrt[((Subscript[A, 1])^2)+((Subscript[A, 2])^2)])cos((Subscript[φ, 1])+3(Subscript[φ, 2]))+(Subscript[H, 5])+...

and we have to compare the values of (Subscript[c, 20])+3(Subscript[c, 11])+9(Subscript[c, 02]) and 3(Sqrt[3(((Subscript[A, 1])^2)+((Subscript[A, 2])^2))]) and apply the Markeev theorem.

The results of our calculations are shown on the picture below.

Conclusion: as the condition |C|=|(Subscript[c, 20])+3(Subscript[c, 11])+9(Subscript[c, 02])|>3(Sqrt[3])A is fulfilled for (Subscript[μ, 1]),(Subscript[μ, 2])  belonging to the fourth order resonance curve between the points A and B, for such values of parameters, according to Markeev's theorem, the equilibrium solution is stable in the Liapunov sense if the fourth order resonance takes place in the system.

The case f=(Subscript[c, 20])((Subscript[σ, 2])^2)+(Subscript[c, 11])(Subscript[σ, 1])(Subscript[σ, 2])+(Subscript[c, 02])((Subscript[σ, 1])^2)=0

There is a curve in the plane (Subscript[μ, 1])(Subscript[Oμ, 2]) where f=0 and conditions of Arnold-Moser theorem are not satisfied. In this case we have to normalize the Hamiltonian up to the sixth order and in the absence of resonances up to the sixth order it takes the form

H=(Subscript[r, 1]) (Subscript[σ, 1])-(Subscript[r, 2]) (Subscript[σ, 2])+(Subscript[c, 20]) (Subscript[r, 1]^2)+(Subscript[c, 02]) (Subscript[r, 2]^2)+(Subscript[c, 11]) (Subscript[r, 1]) (Subscript[r, 2])+(Subscript[c, 30])(Subscript[r, 1]^3)+(Subscript[c, 21])(Subscript[r, 1]^2)(Subscript[r, 2])+(Subscript[c, 12])(Subscript[r, 1])(Subscript[r, 2]^2)+(Subscript[c, 03])(Subscript[r, 2]^3)+(Subscript[H, 7])+...

Our calculations show that the condition

d=(Subscript[c, 30])(Subscript[σ, 2]^3)+(Subscript[c, 21])(Subscript[σ, 2]^2)(Subscript[σ, 1])+(Subscript[c, 12])(Subscript[σ, 2])(Subscript[σ, 1]^2)+(Subscript[c, 03])(Subscript[σ, 1]^3)≠0

is fulfilled for all points where f=0 and we can apply Markeev's theorem.

There are four points of intersection of the curve f=0 and the resonance curves (Subscript[σ, 1])=4(Subscript[σ, 2]), (Subscript[σ, 1])=5(Subscript[σ, 2]) and a question of stability is open for the corresponding values of parameters (Subscript[μ, 1]), (Subscript[μ, 2]).

Conclusion

Summarizing all the results, one can make the following conclusion.

The equilibrium positions (Subscript[S, 5]), (Subscript[S, 6]) in the circular planar restricted problem of four bodies are stable in Liapunov's sense for all values of  the parameters (Subscript[μ, 1]), (Subscript[μ, 2]) from the domain of their linear stability except for 
a) the curve (Subscript[σ, 1])=2(Subscript[σ, 2]), where the third order resonance takes place;
b) the curve (Subscript[σ, 1])=3(Subscript[σ, 2]), where the fourth order resonance takes place, excluding a part of the resonance curve between the points (Subscript[C, 1]) and (Subscript[C, 2]);
c) four points in the plane (Subscript[μ, 1])(Subscript[Oμ, 2]), where the curve f=0 crosses the resonance curves (Subscript[σ, 1])=4(Subscript[σ, 2]) and (Subscript[σ, 1])=5(Subscript[σ, 2]), respectively.

The stability problem at the last four points is still open.

Thank you very much for your attention!!!