Constructing the secular Hamiltonian

Evaluating integrals of multi-pole terms

Here I consider the secular evolution of a test particle interior to a (possibly eccentric) perturber, p. The Hamiltonian governing secular interactions is given by the “double-averaged” potential:
	(Subscript[ℋ, sec])= -(1/(4(π^2)))((∫_0)^(2π))dl((∫_0)^(2π))d (Subscript[l, p])((G (Subscript[m, p]))/(|r-(Subscript[r, p])|))
	= -((G (Subscript[m, p]))/(4(π^2)(Subscript[a, p])))((∫_0)^(2π))dl((∫_0)^(2π))d (Subscript[l, p])  ((∑_(j=2))^∞) (α^j) ((((r/a)))^j)(((((Subscript[a, p])/(Subscript[r, p]))))^(j+1))(Subscript[P, j])[cos(Φ)]
where the second line gives the multi-pole expansion of  the potential in Legendre polynomials. 

I will use a coordinate system in which the x-y plane is the orbital plane of the perturber and the x direction is towards the pericenter of the perturber. In this coordinate system
	cos[Φ] = (1/(1-(Subscript[e, p])cos[(Subscript[u, p])]))(cos[(Subscript[u, p])]-(Subscript[e, p]),(Sqrt[1-((Subscript[e, p])^2)])sin[(Subscript[u, p])],0)·(X,Y,Z)
	 =(((1-(Subscript[e, p])cos[(Subscript[u, p])]))^(-1))[(cos[(Subscript[u, p])]-(Subscript[e, p]))X + (Sqrt[1-((Subscript[e, p])^2)])sin[(Subscript[u, p])]Y]
where (Subscript[u, p]) is the perturber’s eccentric anomaly and (X,Y,Z) are the coordinates of the test particle’s position unit vector (explicit expressions to be given below.)

Integral of (Subscript[l, p])

Integration over (Subscript[l, p]) can be done as follow:
	- First we will use d(Subscript[l, p])=(1-(Subscript[e, p])cos[(Subscript[u, p])])d(Subscript[u, p])to do the integral over eccentric anomaly

	- Note that ((Subscript[r, p])/(Subscript[a, p]))=1-(Subscript[e, p])cos[(Subscript[u, p])]

	- From the expression for cos[Φ] it is clear that the integrand is a series of terms of the form:
		(((1-e cos[(Subscript[u, p])]))^(-n))(cos^(Subscript[k, c]))[(Subscript[u, p])](sin^(Subscript[k, s]))[(Subscript[u, p])]

	- If (Subscript[k, s]) is odd then the integrand term is an odd function and its integral=0.  

	- If (Subscript[k, s]) is even, make the replacement ((sin[(Subscript[u, p])])^2)=1-((cos[(Subscript[u, p])])^2). Thus, all non-zero integrand terms can be expressed in the form:
		 (((1-e cos[(Subscript[u, p])]))^(-n))(cos^k)[(Subscript[u, p])]

Defining 
	 (Subscript[I, n,k])(a,e)=(1/(2π))((∫_0)^(2π))du  (((a-e cos[u]))^(-n))((cos[u])^k)
the recursion relations:
	 (Subscript[I, 1,0])(a,e)=(1/(Sqrt[(a^2)-(e^2)]))
(Subscript[I, n,0])(a,e)=-(1/(n-1))(d/(da))(Subscript[I, n-1,0])(a,e)
(Subscript[I, n,k])(a,e)=(1/(n-1))(d/(de))(Subscript[I, n-1,0])(a,e)
can be used to compute all integral terms by evaluating (Subscript[I, n,k])(1,(Subscript[e, p])).

int[1,0]:=(1/(Sqrt[(a^2)-(ep^2)]));
int[n_Integer,0]:=D[-(1/(n-1))int[n-1,0],a]; 
int[n_Integer,k_Integer]:=(1/(n-1))D[int[n-1,k-1],ep]/;k<n

Below 
	s=sin[(Subscript[u, p])]
 c = cos[(Subscript[u, p])]  
d = 1-(Subscript[e, p])cos[(Subscript[u, p])]

potentialTerm[j_]:=ExpandAll[(d^(-j))LegendreP[j,(((c-ep))/d)X+(Sqrt[1-(ep^2)])(s/d)Y]];

integrationRule={

Times[x__,Power[d,n_],Power[c,k_]]:> Times[x,int[-n,k]/.a->1],
Times[x__,Power[d,n_],c]:> Times[x,int[-n,1]/.a->1],
Times[x__,Power[d,n_]]:> Times[x,int[-n,0]/.a->1]
};

doIntegralOverlp[terms_]:=Simplify[ExpandAll[terms/.(s^2)->1-(c^2)/.s->0]/.integrationRule]

doIntegralOverlp[potentialTerm[2]]

Integral of l

Integration over the particle’s mean longitude is done as follows:
- First, the coordinates of the particle are given by (Subscript[R, z])(h)(Subscript[R, x])(i)(Subscript[R, z])(g)·(((cos[u]-e)/(1-e cos[u])),(((Sqrt[1-(e^2)])sin[u])/(1-e cos[u])),0) where u is the particle’s eccentric anomaly.
- Replacing (r/a)=1-e cos[u] and dl = (1-e cos[u])du, the integrand involves a series of monomials  (((1-e cos[u]))^(j+1))(X^n)(Y^m).
- The integral can be re-written as a series of monomials in sin[u] and cos[u]. 
- Integration over u is then straightforward (Its straightforward to  to verify that no factors  of (1-e cos[u]) will be present in the denominator of the integrand terms, unlike the integral over (Subscript[l, p])

particleCoodinatesRule=Thread[{X,Y,Z}->EulerMatrix[{h,ArcCos[θ],g},{3,1,3}].{((c-e)/(1-e*c)),(Sqrt[1-(e^2)])(s/(1-e*c)),0}];
sincosToComplex = {c->((1/2)(z+(1/z))),s->((1/(2*I))(z-(1/z)))};

integralOverU[expression_]:=Coefficient[ExpandAll[expression/.sincosToComplex],z,0]

multipoleTerm[j_]:=(α^j)Simplify[

integralOverU[(((1-e*c))^(j+1))doIntegralOverlp[potentialTerm[j]]/.particleCoodinatesRule]

]

The Hamiltonian to octupole order

quadrupole=multipoleTerm[2];

octupole=TrigReduce[multipoleTerm[3]];

Compare the expressions below with equations 8-11 in Lithwick+Naoz (2011):

ExpandAll[((8(((1-(ep^2)))^(3/2)))/(3(α^2)))quadrupole]
Simplify[Coefficient[((8(((1-(ep^2)))^(5/2)))/(3(α^3)ep))*(16/5)octupole,{Cos[g+h],Cos[g-h],Cos[3g-h],Cos[3g+h]}]]

Replace orbital elements with dynamical variables, J and Jz

(* Add constant term '(1/3)' to match energy levels listed in LN11, doesn't affect dynamics *)
Fqu = ExpandAll[(((8(((1-(ep^2)))^(3/2)))/(3(α^2)))quadrupole+(1/3))/.θ->((Jz[t])/(Sqrt[1-(e^2)]))/.(e^2)->1-((J[t])^2)/.g->g[t]/.h->h[t]];
Foc =TrigReduce[ Simplify[(((8(((1-(ep^2)))^(5/2)))/(3(α^3)ep))octupole)/.θ->((Jz[t])/(Sqrt[1-(e^2)]))/.e->(Sqrt[1-((J[t])^2)])/.g->g[t]/.h->h[t],Assumptions->J[t]>0]];

Symplectic matrix:

Clear[Jmatrix];
Jmatrix[n_Integer] := ArrayFlatten[{
    {Array[0 &, {n, n}], IdentityMatrix[n]},
    {-IdentityMatrix[n], Array[0 &, {n, n}]}
    }];

Integrate equations of motion

vars={g[t],h[t],J[t],Jz[t]};

With z := (q,p) the equations of motion are given by: ((dz)/(dt)) = J(Subscript[∇, z])ℋ

equationsOfMotion = Thread[D[vars,t]==Jmatrix[2].Grad[Fqu,vars]];
initalCondtions=Thread[(vars/.t->0)=={g0,h0,(Sqrt[1-(e0^2)]),θ0(Sqrt[1-(e0^2)])}];

Specify  initial conditions (setting (Subscript[J, z])=0.2 to match Figure 1 of LN11):

initsrule={g0->Pi/2,h0->0,e0->0.5,θ0->(0.2/(Sqrt[1-(e0^2)]))};

Integrate equations of motion

tStop = 1;
soln=NDSolve[Evaluate[Join[equationsOfMotion,initalCondtions]//.initsrule],(vars/.x_[t]:>x),{t,0,tStop}];

Expression for use in contour plot

Hexprn=Fqu/.{g[t]->g,Jz[t]->0.2,J[t]->0.2/θ};

Show[
(* Plot result of integration *)
ParametricPlot[Evaluate[{g[t],((Jz[t])/(J[t]))}/.soln],{t,0,1},AspectRatio->1,Frame->True,PlotStyle->Directive[Red,Thickness[0.01]]],

(* Contour plot showing level curves of the Hamiltonian *)
ContourPlot[Hexprn,{g,0,2Pi},{θ,0.2,1},Contours->10,ContourShading->None],

(* Plot options*)
PlotRange->{{0,2Pi},{0.0,1}},
Axes->False,
FrameLabel->{ω,θ},
FrameStyle->Directive[16,Black]
]

surface of section with octupole term included

ϵ=0.01;
ℋ=Fqu+ϵ Foc;
equationsOfMotion = Thread[D[vars,t]==Jmatrix[2].Grad[ℋ,vars]];
initalCondtions=Thread[(vars/.t->0)=={g0,h0,(Sqrt[1-(e0^2)]),θ0(Sqrt[1-(e0^2)])}];

replaceRule=(initalCondtions/.Equal->Rule/.x_[0]:>x[t]);

ruleΩPi={g0->0,h0->Pi,e0->e,θ0->θ};
ruleΩ0={g0->0,h0->0,e0->e,θ0->θ};
exprnΩPi=(Fqu)+ϵ Foc/.replaceRule/.ruleΩPi;
exprnΩ0=(Fqu)+ϵ Foc/.replaceRule/.ruleΩ0;

ContourPlot[exprnΩPi,{e,0,0.5},{θ,-1,1},Contours->{0.01,0.04,0.16,0.25,0.36},ContourShading->None,FrameLabel->{e,θ,"Constant Energy Contours"},FrameStyle->Directive[Black,18]]

Section for single trajectory

Integrate equations of motion and record points

Elevel = 0.16

initsrule=ruleΩPi/.First[FindInstance[exprnΩPi==Elevel&&e>0&&θ==0.3,{e,θ}]];

{soln,pts}=Reap@NDSolve[Evaluate[
Join[
equationsOfMotion,
initalCondtions,

{
(* Record ascending node (h) and cosi (Jz/J) every time the argument of perihelion passes g=π *)
WhenEvent[Mod[-g[t],2Pi]==Pi,Sow[{Mod[h[t],2Pi],((Jz[t])/(J[t]))}]]
}
]//.initsrule],
{g,Jz[t],J[t]},{t,0,160}
];

ListPlot[pts,PlotRange->{{0,2Pi},{-1,1}},Frame->True]

GraphicsColumn@{
LogPlot[1-(Sqrt[1-((J[t])^2)])/.soln,{t,0,80},Frame->True,FrameLabel->{t,1-e}],
Plot[((Jz[t])/(J[t]))/.soln,{t,0,80},Frame->True,FrameLabel->{t,θ}],
Plot[Mod[g[t]/.soln,2Pi],{t,0,80},Frame->True,FrameLabel->{t,ω}]
}

Show multiple trajectories

Elevel = 0.16;

allpts=Table[
initsrule=ruleΩPi/.First[FindInstance[exprnΩPi==Elevel&&e>0&&θ==θ0,{e,θ}]];
Last@Last@Reap@NDSolve[Evaluate[
Join[
equationsOfMotion,
initalCondtions,

{
(* Record ascending node (h) and cosi (Jz/J) every time the argument of perihelion passes g=π *)
WhenEvent[Mod[-g[t],2Pi]==Pi,Sow[{Mod[h[t],2Pi],((Jz[t])/(J[t]))}]]
}
]//.initsrule],
{},{t,0,160}
]
,
{θ0,Most@Rest@Subdivide[0,(Sqrt[Elevel]),8]}];

ListPlot[allpts,
Frame->True,FrameLabel->{Ω,θ},
FrameStyle->Directive[Black,18]]

TOI-216

replaceRule=(initalCondtions/.Equal->Rule/.x_[0]:>x[t]);

surface of section with octupole term included

Mstar = 1.1;
Rstar=1.1((Quantity["SolarRadius"])/(Quantity["AstronomicalUnit"]));
Pc=6.2679;
Pb=2093;
eb=0.64;

timeScale = (10*1*^-3)*((2Pi)/Pb)

ab = (Mstar^(1/3))((((Pb/365.25)))^(2/3));
ac = (Mstar^(1/3))((((Pc/365.25)))^(2/3));

ϵ=(((((Pc/Pb)))^(2/3))/(1-(eb^2)))eb;

ℋ=Fqu+ϵ*Foc;

equationsOfMotion = Thread[D[vars,t]==Jmatrix[2].Grad[ℋ,vars]];
initalCondtions=Thread[(vars/.t->0)=={g0,h0,(Sqrt[1-(e0^2)]),θ0(Sqrt[1-(e0^2)])}];

ruleΩPi={g0->0,h0->Pi,e0->e,θ0->θ};
ruleΩ0={g0->0,h0->0,e0->e,θ0->θ};
exprnΩPi=(Fqu)+ϵ Foc/.replaceRule/.ruleΩPi;
exprnΩ0=(Fqu)+ϵ Foc/.replaceRule/.ruleΩ0;

ContourPlot[exprnΩPi,{e,0,0.5},{θ,-1,1},Contours->{0.01,0.04,0.16,0.25,0.36},ContourShading->None,FrameLabel->{e,θ,"Constant Energy Contours"},FrameStyle->Directive[Black,18]]

Fqu/.{g[t]->0,h[t]->Pi,Jz[t]->θ0(Sqrt[1-(e0^2)]),J[t]->(Sqrt[1-(e0^2)])}/.e0->0

Section for single trajectory

Integrate equations of motion and record points

tstop=2 160;
initsrule={g0->Pi/2,h0->0,e0->0.001,θ0->Cos[50*Degree]};
{soln,pts}=Reap@NDSolve[Evaluate[
Join[
equationsOfMotion,
initalCondtions,

{
(* Record ascending node (h) and cosi (Jz/J) every time the argument of perihelion passes g=π *)
WhenEvent[Mod[-g[t],2Pi]==Pi,Sow[{Mod[h[t],2Pi],((Jz[t])/(J[t]))}]]
}
]//.initsrule],
{g,Jz[t],J[t]},{t,0,tstop}
];
GraphicsColumn@{
LogPlot[(ac/Rstar)(1-(Sqrt[1-((J[t])^2)]))/.soln,{t,0,tstop},Frame->True,FrameLabel->{t,a(1-e)/(SubStar[R])},PlotRange->{3,20}],
Plot[((Jz[t])/(J[t]))/.soln,{t,0,tstop},Frame->True,FrameLabel->{t,θ}],
Plot[Mod[g[t]/.soln,2Pi],{t,0,tstop},Frame->True,FrameLabel->{t,ω}]
}

ListPlot[pts,PlotRange->{{0,2Pi},{-1,1}},Frame->True]

Show multiple trajectories

Elevel = 0.16;

allpts=Table[
initsrule=ruleΩPi/.First[FindInstance[exprnΩPi==Elevel&&e>0&&θ==θ0,{e,θ}]];
Last@Last@Reap@NDSolve[Evaluate[
Join[
equationsOfMotion,
initalCondtions,

{
(* Record ascending node (h) and cosi (Jz/J) every time the argument of perihelion passes g=π *)
WhenEvent[Mod[-g[t],2Pi]==Pi,Sow[{Mod[h[t],2Pi],((Jz[t])/(J[t]))}]]
}
]//.initsrule],
{},{t,0,160}
]
,
{θ0,Most@Rest@Subdivide[0,(Sqrt[Elevel]),8]}];

ListPlot[allpts,
Frame->True,FrameLabel->{Ω,θ},
FrameStyle->Directive[Black,18]]

initalCondtions=Thread[(vars/.t->0)=={g0,h0,(Sqrt[1-(e0^2)]),θ0(Sqrt[1-(e0^2)])}];
initsrule={g0->Pi/2,h0->0,e0->0.5,θ0->(0.2/(Sqrt[1-(e0^2)]))};

GraphicsColumn@{
LogPlot[1-(Sqrt[1-((J[t])^2)])/.soln,{t,0,80},Frame->True,FrameLabel->{t,1-e}],
Plot[((Jz[t])/(J[t]))/.soln,{t,0,80},Frame->True,FrameLabel->{t,θ}],
Plot[Mod[g[t]/.soln,2Pi],{t,0,80},Frame->True,FrameLabel->{t,ω}]
}