Collinear Three-Body Problem

This notebook explores the collinear three-body problem using McGehee coordinates.  After executing the first section you can just go on to look at examples.

Coordinate changes and solving the ODEs

This section contains functions relating Cartesian, Jacobi and Blown-Up coordinates as well as functions to calculate an orbit of the collinear 3bp, find and set energy, etc.

The positions of the bodies will be q1 < q3 < q2 along the line.  If the center of mass is at the origin, we can replace these three coordinates by two Jacobi variables.  The velocities transform in the same way.

CartesianToJacobi[{q1_,q2_,q3_},{m1_,m2_,m3_}]:=With[{a1=m1/(m1+m2),a2=m2/(m1+m2)},{q2-q1,q3-(a1*q1+a2*q2)}]

CartesianToJacobi[{q1_,q2_,q3_,v1_,v2_,v3_},{m1_,m2_,m3_}]:=With[{a1=m1/(m1+m2),a2=m2/(m1+m2)},{q2-q1,q3-(a1*q1+a2*q2),v2-v1,v3-(a1*v1+a2*v2)}]

 JacobiToCartesian[{x1_,x2_},{m1_,m2_,m3_}]:=With[{a1=m1/(m1+m2),a2=m2/(m1+m2)},{-a2*x1-m3*x2/(m1+m2+m3),a1*x1-m3*x2/(m1+m2+m3),(m1+m2)*x2/(m1+m2+m3)}]

JacobiToCartesian[{x1_,x2_,u1_,u2_},{m1_,m2_,m3_}]:=With[{a1=m1/(m1+m2),a2=m2/(m1+m2)},{-a2*x1-m3*x2/(m1+m2+m3),a1*x1-m3*x2/(m1+m2+m3),(m1+m2)*x2/(m1+m2+m3),-a2*u1-m3*u2/(m1+m2+m3),a1*u1-m3*u2/(m1+m2+m3),(m1+m2)*u2/(m1+m2+m3)}]

Next introduce blown-up coordinates (r,theta,nu,gamma).

JacobiToBlowUp[{x1_,x2_},{m1_,m2_,m3_}]:=With[{mu1=m1*m2/(m1+m2),mu2=(m1+m2)*m3/(m1+m2+m3)},{Sqrt[mu1*x1^2+mu2*x2^2],ArcTan[Sqrt[mu1]*x1,Sqrt[mu2]*x2]}]

JacobiToBlowUp[{x1_,x2_,u1_,u2_},{m1_,m2_,m3_}]:=Module[{a1=m1/(m1+m2),a2=m2/(m1+m2),mu1=m1*m2/(m1+m2),mu2=(m1+m2)*m3/(m1+m2+m3),r,nu,A1,A2,phi},
A1=Sqrt[(m1+m3)/(m1*m3)];
A2=Sqrt[(m2+m3)/(m2*m3)];
phi = (x2+a2*x1)(x2-a1*x1)/(A1*A2);
r=Sqrt[mu1*x1^2+mu2*x2^2];
nu=(mu1*x1*u1+mu2*x2*u2)/Sqrt[r];{r,ArcTan[Sqrt[mu1]*x1,Sqrt[mu2]*x2],nu, phi*Sqrt[mu1*mu2/r]*(x1*u2-x2*u1)}]

BlowUpToJacobi[{r_,theta_},{m1_,m2_,m3_}]:=With[{mu1=m1*m2/(m1+m2),mu2=(m1+m2)*m3/(m1+m2+m3)},{r*Cos[theta]/Sqrt[mu1],r*Sin[theta]/Sqrt[mu2]}]

BlowUpToJacobi[{r_,theta_,nu_,gamma_},{m1_,m2_,m3_}]:=Module[{a1=m1/(m1+m2),a2=m2/(m1+m2),mu1=m1*m2/(m1+m2),mu2=(m1+m2)*m3/(m1+m2+m3),x1,x2,A1,A2,phi},
A1=Sqrt[(m1+m3)/(m1*m3)];
A2=Sqrt[(m2+m3)/(m2*m3)];
{x1,x2}={r*Cos[theta]/Sqrt[mu1],r*Sin[theta]/Sqrt[mu2]};
phi = (x2+a2*x1)(x2-a1*x1)/(A1*A2);
{x1,x2,(x1*nu-Sqrt[mu2/mu1]*x2*gamma/phi)/(r^(3/2)),(x2*nu+Sqrt[mu1/mu2]*x1*gamma/phi)/(r^(3/2))}]

CartesianToBlowUp[{q1_,q2_,q3_},{m1_,m2_,m3_}]:=JacobiToBlowUp[CartesianToJacobi[{q1,q2,q3},{m1,m2,m3}],{m1,m2,m3}]

CartesianToBlowUp[{q1_,q2_,q3_,v1_,v2_,v3_},{m1_,m2_,m3_}]:=JacobiToBlowUp[CartesianToJacobi[{q1,q2,q3,v1,v2,v3},{m1,m2,m3}],{m1,m2,m3}]

BlowUpToCartesian[{r_,theta_},{m1_,m2_,m3_}]:=JacobiToCartesian[BlowUpToJacobi[{r,theta},{m1,m2,m3}],{m1,m2,m3}]

BlowUpToCartesian[{r_,theta_,nu_,gamma_},{m1_,m2_,m3_}]:=JacobiToCartesian[BlowUpToJacobi[{r,theta,nu,gamma},{m1,m2,m3}],{m1,m2,m3}]

The  regularized  “potential”, W(theta), its derivative and the time rescaling factor phi

V1=m1*m2/r12+m1*m3/r13+m2*m3/r23
V=V1/.{r12->Cos[theta]/Sqrt[mu1],r13->Sin[theta]/Sqrt[mu2]+a2*Cos[theta]/Sqrt[mu1],r23->-Sin[theta]/Sqrt[mu2]+a1*Cos[theta]/Sqrt[mu1]}/.{a1->m1/(m1+m2),a2->m2/(m1+m2),mu1->m1*m2/(m1+m2),mu2->(m1+m2)*m3/(m1+m2+m3)}
Vtheta = D[V,theta];


W=W1/.{r12->Cos[theta]/Sqrt[mu1],r13->Sin[theta]/Sqrt[mu2]+a2*Cos[theta]/Sqrt[mu1],r23->-Sin[theta]/Sqrt[mu2]+a1*Cos[theta]/Sqrt[mu1]}/.{a1->m1/(m1+m2),a2->m2/(m1+m2),mu1->m1*m2/(m1+m2),mu2->(m1+m2)*m3/(m1+m2+m3)}/.{A1->Sqrt[(m1+m3)/(m1*m3)],
A2->Sqrt[(m2+m3)/(m2*m3)]};
Wtheta = D[W,theta];

W1=m1*m2*r13*r23/(A1*A2*r12)+m1*m3*r23/(A1*A2)+m2*m3*r13/(A1*A2)
W=W1/.{r12->Cos[theta]/Sqrt[mu1],r13->Sin[theta]/Sqrt[mu2]+a2*Cos[theta]/Sqrt[mu1],r23->-Sin[theta]/Sqrt[mu2]+a1*Cos[theta]/Sqrt[mu1]}/.{a1->m1/(m1+m2),a2->m2/(m1+m2),mu1->m1*m2/(m1+m2),mu2->(m1+m2)*m3/(m1+m2+m3)}/.{A1->Sqrt[(m1+m3)/(m1*m3)],
A2->Sqrt[(m2+m3)/(m2*m3)]};
Wtheta = D[W,theta];
phi1=r13*r23/(A1*A2)
phi = phi1/.{r12->Cos[theta]/Sqrt[mu1],r13->Sin[theta]/Sqrt[mu2]+a2*Cos[theta]/Sqrt[mu1],r23->-Sin[theta]/Sqrt[mu2]+a1*Cos[theta]/Sqrt[mu1]}/.{a1->m1/(m1+m2),a2->m2/(m1+m2),mu1->m1*m2/(m1+m2),mu2->(m1+m2)*m3/(m1+m2+m3)}/.{A1->Sqrt[(m1+m3)/(m1*m3)],
A2->Sqrt[(m2+m3)/(m2*m3)]};
phitheta=D[phi,theta];

Calculate  the indpoints of the theta interval.   The endpoints correspond to binary collisions.

ThetaBounds[{m10_,m20_,m30_}]:=Module[{a1=m10/(m10+m20),a2=m20/(m10+m20),mu1=m10*m20/(m10+m20),mu2=(m10+m20)*m30/(m10+m20+m30)},
{ArcTan[-a2*Sqrt[mu2/mu1]],ArcTan[a1*Sqrt[mu2/mu1]]}]

Calculate the energy, h, in McGehee coordinates.  Won't work at the collision points.  The regularized energy is the converted energy function  which should be zero.

Energy[{r0_,theta0_,nu0_,gamma0_,___},{m10_,m20_,m30_}]:=Module[{phi,a1,a2,mu1,mu2},
a1=m10/(m10+m20);a2=m20/(m10+m20);mu1=m10*m20/(m10+m20);mu2=(m10+m20)*m30/(m10+m20+m30);
phi=(Sin[theta0]/Sqrt[mu2]+a2*Cos[theta0]/Sqrt[mu1])(-Sin[theta0]/Sqrt[mu2]+a1*Cos[theta0]/Sqrt[mu1])/(Sqrt[(m10+m30)/(m10*m30)]*Sqrt[(m20+m30)/(m20*m30)]);
(phi^2*nu0^2/2+gamma0^2/2-phi*W)/(r0*phi^2)/.{theta->theta0,m1->m10,m2->m20,m3->m30}]

(* The following energy function should always be zero *)
RegularizedEnergy[{r0_,theta0_,nu0_,gamma0_,___},{m10_,m20_,m30_},h0_]:=Module[{phi,a1,a2,mu1,mu2},
a1=m10/(m10+m20);a2=m20/(m10+m20);mu1=m10*m20/(m10+m20);mu2=(m10+m20)*m30/(m10+m20+m30);
phi=(Sin[theta0]/Sqrt[mu2]+a2*Cos[theta0]/Sqrt[mu1])(-Sin[theta0]/Sqrt[mu2]+a1*Cos[theta0]/Sqrt[mu1])/(Sqrt[(m10+m30)/(m10*m30)]*Sqrt[(m20+m30)/(m20*m30)]);
(phi^2*nu0^2/2+gamma0^2/2-phi*W)-r0*phi^2(h0)/.{theta->theta0,m1->m10,m2->m20,m3->m30}]

Regularized ODE

rprime = phi*nu*r;
thetaprime=gamma;
nuprime = W-phi*nu^2/2+2*r*h*phi;
gammaprime = phi*Wtheta-phi*nu*gamma/2 +phitheta*(2*r*h*phi+W-phi*nu^2);
timeprime=phi;

collinearvf={rprime,thetaprime, nuprime, gammaprime, timeprime};

A function to solve the regularized ODE  starting from a given initial condition. First it calculates the energy so it can substitute that into the regularized ODE.

vf=collinearvf/.{r->r[t],theta->theta[t],nu->nu[t],gamma->gamma[t]};
MakeRegularizedOrbit[{r0_,theta0_,nu0_,gamma0_},{m10_,m20_,m30_},tmax_]:=
Module[{vf2,h0},
h0=If[r0>0,Energy[{r0,theta0,nu0,gamma0},{m10,m20,m30}],-1];
vf2=vf/.{m1->m10,m2->m20,m3->m30,h->h0};NDSolveValue[{r'[t]==vf2[[1]],theta'[t]==vf2[[2]],nu'[t]==vf2[[3]],gamma'[t]==vf2[[4]], time'[t]==vf2[[5]],r[0]==r0,theta[0]==theta0, nu[0]==nu0,gamma[0]==gamma0,time[0]==0},{r[t],theta[t],nu[t],gamma[t], time[t]},{t,0,tmax}]]


MakeRegularizedOrbitFromCartesianIC[{q10_,q20_,q30_,v10_,v20_,v30_},{m10_,m20_,m30_},tmax_]:=
Module[{vf2,h0,r0,theta0,nu0,gamma0},
{r0,theta0,nu0,gamma0} =CartesianToBlowUp[{q10,q20,q30,v10,v20,v30},{m10,m20,m30}];
h0=If[r0>0,Energy[{r0,theta0,nu0,gamma0},{m10,m20,m30}],-1];
vf2=vf/.{m1->m10,m2->m20,m3->m30,h->h0};NDSolveValue[{r'[t]==vf2[[1]],theta'[t]==vf2[[2]],nu'[t]==vf2[[3]],gamma'[t]==vf2[[4]], time'[t]==vf2[[5]],r[0]==r0,theta[0]==theta0, nu[0]==nu0,gamma[0]==gamma0,time[0]==0},{r[t],theta[t],nu[t],gamma[t], time[t]},{t,0,tmax}]]

Set the energy starting from some initial condition in blown-up coordinates.  This can be done by resetting one or the other the coordinates.  The generic version will be done by resetting r unless r=0 in which case, it resets gamma or nu.

Clear[SetEnergyByNu]

SetEnergy[{r0_,theta0_,nu0_,gamma0_},{m10_,m20_,m30_},h0_]:=Module[{newgammasq,phi,a1,a2,mu1,mu2},
If[r0>0,Return[SetEnergyByR[{r0,theta0,nu0,gamma0},{m10,m20,m30},h0]]];
If[gamma0^2>0.0001,Return[SetEnergyByGamma[{r0,theta0,nu0,gamma0},{m10,m20,m30},h0]]];
SetEnergyByNu[{r0,theta0,nu0,gamma0},{m10,m20,m30},h0]
]

SetEnergyByR[{r0_,theta0_,nu0_,gamma0_},{m10_,m20_,m30_},h0_]:=Module[{newr,phi,a1,a2,mu1,mu2},
a1=m10/(m10+m20);a2=m20/(m10+m20);mu1=m10*m20/(m10+m20);mu2=(m10+m20)*m30/(m10+m20+m30);
phi=(Sin[theta0]/Sqrt[mu2]+a2*Cos[theta0]/Sqrt[mu1])(-Sin[theta0]/Sqrt[mu2]+a1*Cos[theta0]/Sqrt[mu1])/(Sqrt[(m10+m30)/(m10*m30)]*Sqrt[(m20+m30)/(m20*m30)]);
newr=(phi^2*nu0^2/2+gamma0^2/2-phi*W)/(h0*phi^2)/.{theta->theta0,m1->m10,m2->m20,m3->m30};
{newr,theta0,nu0,gamma0}]

SetEnergyByGamma[{r0_,theta0_,nu0_,gamma0_},{m10_,m20_,m30_},h0_]:=Module[{newgammasq,phi,a1,a2,mu1,mu2},
a1=m10/(m10+m20);a2=m20/(m10+m20);mu1=m10*m20/(m10+m20);mu2=(m10+m20)*m30/(m10+m20+m30);
phi=(Sin[theta0]/Sqrt[mu2]+a2*Cos[theta0]/Sqrt[mu1])(-Sin[theta0]/Sqrt[mu2]+a1*Cos[theta0]/Sqrt[mu1])/(Sqrt[(m10+m30)/(m10*m30)]*Sqrt[(m20+m30)/(m20*m30)]);
newgammasq=2*r0*(h0)*phi^2+2*phi*W-phi^2*nu0^2/.{theta->theta0,m1->m10,m2->m20,m3->m30};
{r0,theta0,nu0,Sign[gamma0]*Sqrt[newgammasq]}]

SetEnergyByNu[{r0_,theta0_,nu0_,gamma0_},{m10_,m20_,m30_},h0_]:=Module[{newnusq,phi,a1,a2,mu1,mu2},
a1=m10/(m10+m20);a2=m20/(m10+m20);mu1=m10*m20/(m10+m20);mu2=(m10+m20)*m30/(m10+m20+m30);
phi=(Sin[theta0]/Sqrt[mu2]+a2*Cos[theta0]/Sqrt[mu1])(-Sin[theta0]/Sqrt[mu2]+a1*Cos[theta0]/Sqrt[mu1])/(Sqrt[(m10+m30)/(m10*m30)]*Sqrt[(m20+m30)/(m20*m30)]);
newnusq=(2*r0*(h0)*phi^2+2*phi*W-gamma0^2)/phi^2/.{theta->theta0,m1->m10,m2->m20,m3->m30};
{r0,theta0,Sign[nu0]*Sqrt[newnusq],gamma0}]

Examples and plots

The simplest way to compute an orbit is just to enter initial values for the masses and the McGehee coordinates {r,theta,nu, omega}.   You can calculate the admissible interval {θ1,θ2} and the energy, and calculate an orbit for a given length of time.

m10=1;m20=2;m30=3;
r0=1;
theta0=0;
nu0=0;
gamma0=0;
tmax0=10;
tb0=ThetaBounds[{m10,m20,m30}]
h0=Energy[{r0,theta0,nu0,gamma0},{1,2,3}]//N
orb=MakeRegularizedOrbit[{r0,theta0,nu0,gamma0},{m10,m20,m30},tmax0];

Here is a plot of  the shape θ and size r with the bounds for theta (the double collisions).

bounds=Graphics[{Line[{{tb0[[1]],0},{tb0[[1]],10}}],Line[{{tb0[[2]],0},{tb0[[2]],10}}]}];
Show[ParametricPlot[orb[[{2,1}]],{t,0,10},AxesLabel->{θ,r}],bounds]

To see what the three  bodies are doing, it’s nice to have the following functions.  The first one works with Cartesian coordinates, the second with McGehee coordinates, and the last just plots the shape.

MakeCollinearPlot[{q1_,q2_,q3_},{m1_,m2_,m3_},width_:3]:=Module[{p1,p2,p3,l,ps,},
ps= 0.03;
p1=Graphics[{Blue,PointSize[ps],Point[{q1,0}]}];
p2=Graphics[{Green,PointSize[ps],Point[{q2,0}]}];
p3=Graphics[{Black,PointSize[ps],Point[{q3,0}]}];
l=Graphics[Line[{{-width,0},{width,0}}]];
Show[p1,p2,p3,l,PlotRange->{{-width,width},{-0.5,0.5}}]]

MakeCollinearPlot[{r_,theta_},{m1_,m2_,m3_},width_]:=Module[{q1,q2,q3},
{q1,q2,q3}=BlowUpToCartesian[{r,theta},{m1,m2,m3}];
MakeCollinearPlot[{q1,q2,q3},{m1,m2,m3},width]]

MakeCollinearPlot[theta_,{m1_,m2_,m3_},width_:1.2]:=Module[{q1,q2,q3},
{q1,q2,q3}=BlowUpToCartesian[{1,theta},{m1,m2,m3}];
MakeCollinearPlot[{q1,q2,q3},{m1,m2,m3},width]]

Here is an  animation of the orbit calculated above.

Manipulate[MakeCollinearPlot[orb[[{1,2}]]/.t->t0,{m10,m20,m30},3],{t0,0,tmax0}]

You can also start with a Cartesian initial condition which will be converted to have center of mass at 0.  Positions should have q1<q3<q2

m10=1;m20=2;m30=2;
q10=1;q20=5;q30=2;
v10=0;v20=0;v30=0;
tmax0=20;
tb0=ThetaBounds[{m10,m20,m30}]//N
orb=MakeRegularizedOrbitFromCartesianIC[{q10,q20,q30,v10,v20,v30},{m10,m20,m30},tmax0];
Manipulate[MakeCollinearPlot[orb[[{1,2}]]/.t->t0,{m10,m20,m30},5],{t0,0,tmax0}]

Another interesting plot is a spacetime plot showing the three positions versus time.

{f1,f2,f3}=BlowUpToCartesian[orb[[{1,2}]],{m10,m20,m30}];
spacetimeplot =Plot[{f1,f2,f3},{t,0,11},AxesLabel->{t,q},PlotStyle->{Blue,Green,Black}];
t0=0;
pts[t0_]:=Graphics[{PointSize[0.02],Blue, Point[{t0,f1/.t->t0}],Green, Point[{t0,f2/.t->t0}],Black, Point[{t0,f3/.t->t0}]}];
pttable=Table[pts[t0],{t0,0,11}];
Show[spacetimeplot,pttable]

Central Configuration and Triple Collision Orbits

The Euler central configuration is the critical point of the unregularized potential V = W/phi

m10=1;m20=2;m30=3;
tb0=ThetaBounds[{m10,m20,m30}];
Plot[W/phi/.{m1->m10,m2->m20,m3->m30},{theta,tb0[[1]],tb0[[2]]}]

The  following solves numerically for the position of the CC and of  the corresponding restpoints on the collision manifold.  There are two restpoints with different signs for nu.

FindCC[{m10_,m20_,m30_}]:=theta/.FindRoot[Vtheta/.{m1->m10,m2->m20,m3->m30},{theta,0}][[1]]

FindEulerRestpoint[{m10_,m20_,m30_},sign_]:=With[{theta0=FindCC[{m10,m20,m30}]},
{0,theta0,sign*Sqrt[2*W/phi]/.{theta->theta0,m1->m10,m2->m20,m3->m30},0}]

FindCC[{1,2,3}]
FindEulerRestpoint[{1,2,3},1]

Here is the collision manifold for these masses projected to {θ,γ,ν} space with the level curves of the Lyapunov function ν and the restpoints.

eulerrpminus = FindEulerRestpoint[{1,2,3},-1];
eulerrpplus = FindEulerRestpoint[{1,2,3},1];
restpoints =Graphics3D[{PointSize[0.03],Point[eulerrpminus[[{2,4,3}]]],Point[eulerrpplus[[{2,4,3}]]]}];
F = phi^2nu^2/2+gamma^2/2-phi*W /.{m1->1,m2->2,m3->3};
collmanplot = ContourPlot3D[F==0,{theta,tb0[[1]],tb0[[2]]},{gamma,-5,5},{nu,-10,10},MeshFunctions->{#3&},ContourStyle->Opacity[0.5]];
collmanplot2=Show[collmanplot,restpoints]

Manifolds of triple collision orbits.

Next we will compute approximations to the stable and unstable manifolds of the  restpoints.  The 2D unstable manifold of the  restpoint  p+ with  ν > 0 gives the orbits which have triple collision in backward time.  The 2D stable manifold of the  restpoint  p- with  ν < 0 gives the orbits which have triple collision in forward time. We will compute the unstable manifold by forward time integration starting from initial conditions close to the unstable  eigenspace.  Here are the restpoints, eigenvalues and eigenvectors for masses 1,2,3.

eulerrpminus = FindEulerRestpoint[{1,2,3},-1]
eulerrpplus = FindEulerRestpoint[{1,2,3},1];
vf0=collinearvf[[{1,2,3,4}]]/.{m1->1,m2->2,m3->3,h->-1}//Simplify;
dvf0 = Map[D[vf0,#]&,{r,theta,nu,gamma}]//Transpose;
vf0/.{r->eulerrpplus[[1]],theta->eulerrpplus[[2]],nu->eulerrpplus[[3]],gamma->eulerrpplus[[4]]};
dvf0=dvf0/.{r->eulerrpplus[[1]],theta->eulerrpplus[[2]],nu->eulerrpplus[[3]],gamma->eulerrpplus[[4]]};
%//MatrixForm
evals=Eigenvalues[dvf0]
evects=Eigenvectors[dvf0]

The relevant  eigenvectors are numbers 2 and 3 corresponding to the positive eigenvalues.  Here is a function to make initial conditions in this subspace.   The parameter called angle is to control the blend of the two eigenvectors.  Then we make a table of these initial conditions.

MakeUnstableIC[angle_]:=SetEnergyByNu[eulerrpplus + 0.01*Cos[angle]*evects[[3]]+0.00001*Sin[angle]*evects[[2]],{1,2,3},-1]

wutable=Table[MakeUnstableIC[angle],{angle,-Pi/2,Pi/2,Pi/12}]

Here is a nice plot of the resulting orbits.

tb0=ThetaBounds[{1,2,3}];
wuorbits=Map[MakeRegularizedOrbit[#,{1,2,3},10]&,wutable];
wuplotplus=ParametricPlot[Map[Part[#,{2,3}]&,wuorbits],{t,0,10},PlotRange->{{tb0[[1]],tb0[[2]]},{-6,10}},AspectRatio->1];
F= phi^2*nu^2/2-W*phi/.{m1->1,m2->2,m3->3};
cp=ContourPlot[F, {theta,tb0[[1]],tb0[[2]]},{nu,-10,10},Contours->{0},ContourShading->{LightBlue,None},FrameLabel->{θ,ν}];
Show[cp,wuplotplus,Graphics[{PointSize[0.03],Point[eulerrpminus[[{2,3}]]],Point[eulerrpplus[[{2,3}]]]}]]

And here is a 3D plot.  Evidently this 2D unstable manifold  is a very complicated surface in the solid region inside the collision manifold.

wuplotplus3D=ParametricPlot3D[Map[Part[#,{2,4,3}]&,wuorbits],{t,0,10},PlotRange->{{tb0[[1]]-0.1,tb0[[2]]+0.1},{-5,5}, {-6,10}},BoxRatios->{10,10,16}];
F = phi^2nu^2/2+gamma^2/2-phi*W /.{m1->1,m2->2,m3->3};
collmanplot3 = ContourPlot3D[F==0,{theta,tb0[[1]],tb0[[2]]},{gamma,-5,5},{nu,4,8},MeshFunctions->{#3&},ContourStyle->Opacity[0.2]];
Show[wuplotplus3D,collmanplot3]

Saddle point branches in the collision manifold and nonregularizability

Next we will compute approximations to the stable and unstable manifolds of the  restpoints.  For the restpoint with ν < 0, the 2D stable manifold gives the forward time triple collision orbits while  the 1D unstable manifold lies entirely in the collision manifold but plays an important  role  for  near-collision orbits.  Here are the restpoints, eigenvalues and eigenvectors for masses 1,2,3.

eulerrpminus = FindEulerRestpoint[{1,2,3},-1]
eulerrpplus = FindEulerRestpoint[{1,2,3},1];
vf0=collinearvf[[{1,2,3,4}]]/.{m1->1,m2->2,m3->3,h->-1}//Simplify;
dvf0 = Map[D[vf0,#]&,{r,theta,nu,gamma}]//Transpose;
vf0/.{r->eulerrpminus[[1]],theta->eulerrpminus[[2]],nu->eulerrpminus[[3]],gamma->eulerrpminus[[4]]}
dvf0=dvf0/.{r->eulerrpminus[[1]],theta->eulerrpminus[[2]],nu->eulerrpminus[[3]],gamma->eulerrpminus[[4]]}
%//MatrixForm
evals=Eigenvalues[dvf0]
evects=Eigenvectors[dvf0]

The last eigenvector should give a good approximation to the branches of the saddle with ν<0.

eulerrpminus+0.01*evects[[4]]
ic = SetEnergyByNu[eulerrpminus+0.01*evects[[1]],{1,2,3},-1]
ic = SetEnergyByNu[eulerrpminus+{0,0.01,0,0.01},{1,2,3},-1]
wuorb1=MakeRegularizedOrbit[ic,{1,2,3},10];
eulerrpminus+0.01*evects[[4]]
ic = SetEnergyByNu[eulerrpminus-0.01*evects[[1]],{1,2,3},-1]
ic = SetEnergyByNu[eulerrpminus-{0,0.01,0,0.01},{1,2,3},-1]
wuorb2=MakeRegularizedOrbit[ic,{1,2,3},10];
Plot[{wuorb1[[2]],wuorb2[[2]]},{t,0,10}]

A 2D plot in the  {θ,ν} plane is shows the behavior nicely.

F= phi^2*nu^2/2-W*phi/.{m1->1,m2->2,m3->3};
cp=ContourPlot[F, {theta,tb0[[1]],tb0[[2]]},{nu,-10,10},Contours->{0},ContourShading->{LightBlue,None},FrameLabel->{θ,ν}];
wuplot  =ParametricPlot[{wuorb1[[{2,3}]],wuorb2[[{2,3}]]},{t,0,10},PlotRange->{tb0,{-8,12}},AspectRatio->1,AxesLabel->{θ,ν},PlotStyle->{Blue,Green}];
Show[cp,wuplot,Graphics[{PointSize[0.03],Point[eulerrpminus[[{2,3}]]],Point[eulerrpplus[[{2,3}]]]}]]

The stable manifolds of the restpoint with ν >0  are obtained by symmetry.

Reflect[{r0_,theta0_,nu0_,gamma0_,___}]:={r0,theta0,-nu0,-gamma0}
wsorb1=Reflect[wuorb1];
wsorb2=Reflect[wuorb2];wsplot  =ParametricPlot[{wsorb1[[{2,3}]],wsorb2[[{2,3}]]},{t,0,10},PlotRange->{tb0,{-8,12}},AspectRatio->1,AxesLabel->{θ,ν},PlotStyle->{Orange,Red}];
Show[cp,wuplot,wsplot,Graphics[{PointSize[0.03],Point[eulerrpminus[[{2,3}]]],Point[eulerrpplus[[{2,3}]]]}]]

One of the interesting questions about the flow on the triple collision  manifold is:  how does the behavior of the branches of the saddle points depend on the choice of masses.  To  explore this, we can  package up some of the methods  used for the special  masses 1,2, 3 above.

FindSaddleManifolds[{m10_,m20_,m30_},tmax_:10]:= Module[{eplus, eminus, wu1, wu2,dv,F,cp,wuplot,tb0,ic,wuorb1,wuorb2},
eplus = FindEulerRestpoint[{m10,m20,m30},1];
eminus = FindEulerRestpoint[{m10,m20,m30},-1];
dv = 2*eplus[[3]];  (* a good size for the plots *)
tb0=ThetaBounds[{m10,m20,m30}];
F= phi^2*nu^2/2-W*phi/.{m1->m10,m2->m20,m3->m30};
cp=ContourPlot[F, {theta,tb0[[1]],tb0[[2]]},{nu,-dv,dv},Contours->{0},ContourShading->{LightBlue,None},FrameLabel->{θ,ν}];
ic = SetEnergyByNu[eminus+{0,0.01,0,0.01},{m10,m20,m30},-1/2];
wuorb1=MakeRegularizedOrbit[ic,{m10,m20,m30},tmax];
ic = SetEnergyByNu[eminus-{0,0.01,0,0.01},{m10,m20,m30},-1/2];
wuorb2=MakeRegularizedOrbit[ic,{m10,m20,m30},tmax];
wuplot  =ParametricPlot[{wuorb1[[{2,3}]],wuorb2[[{2,3}]]},{t,0,tmax},PlotRange->{tb0,{-dv,dv}},AspectRatio->1,AxesLabel->{θ,ν},PlotStyle->{Blue,Green}];
Show[cp,wuplot,Graphics[{PointSize[0.03],Point[eminus[[{2,3}]]],Point[eplus[[{2,3}]]]}]]]

Manipulate[FindSaddleManifolds[{1,2,m30},100],{{m30,1},0,3}]

s01=FindSaddleManifolds[{1,2,0.1},200];
s1=FindSaddleManifolds[{1,2,1},20];
s2=FindSaddleManifolds[{1,2,2},20];
s3=FindSaddleManifolds[{1,2,3},20];

Show[GraphicsGrid[{{s01,s1},{s2,s3}}]]

Two nearby near collision orbits with very different fates, illustrating the nonregularizability of triple collision. Start in the plane ν=0 at and close to the homothetic solution.

EulerCC= FindCC[{1,2,3}];
eps = 0.00001;
h0=Energy[{1,EulerCC,0,0},{1,2,3}]
ic0 = SetEnergy[{1,EulerCC,0,0},{1,2,3},h0];
ic1 = SetEnergy[ic0+{0,eps,0,0},{1,2,3},h0];
ic2 = SetEnergy[ic0-{0,eps,0,0},{1,2,3},h0];
orb0=MakeRegularizedOrbit[ic1,{1,2,3},20];
orb1=MakeRegularizedOrbit[ic1,{1,2,3},20];
orb2=MakeRegularizedOrbit[ic2,{1,2,3},20];
{f1,f2,f3}=BlowUpToCartesian[orb0[[{1,2}]],{1,2,3}];
p0=Plot[{f1,f2,f3},{t,0,20},AxesLabel->{t,q},PlotStyle->{Blue,Green,Black},PlotRange->All];
{f1,f2,f3}=BlowUpToCartesian[orb1[[{1,2}]],{1,2,3}];
p1=Plot[{f1,f2,f3},{t,0,20},AxesLabel->{t,q},PlotStyle->{Blue,Green,Black},PlotRange->All];
{f1,f2,f3}=BlowUpToCartesian[orb2[[{1,2}]],{1,2,3}];
p2=Plot[{f1,f2,f3},{t,0,20},AxesLabel->{t,q},PlotStyle->{Blue,Green,Black},PlotRange->All];
Show[GraphicsGrid[{{p1,p2}}]]

{f1,f2,f3}=BlowUpToCartesian[orb1[[{1,2}]],{1,2,3}];
p1=Plot[{f1,f2,f3},{t,3,8},AxesLabel->{t,q},PlotStyle->{Blue,Green,Black},PlotRange->All];
{f1,f2,f3}=BlowUpToCartesian[orb2[[{1,2}]],{1,2,3}];
p2=Plot[{f1,f2,f3},{t,3,8},AxesLabel->{t,q},PlotStyle->{Blue,Green,Black},PlotRange->All];
Show[GraphicsGrid[{{p1,p2}}]]