(* ::Package:: *)

(* 
This file is associated with "Scattering Amplitudes, the Tail Effect, 
and Conservative Binary Dynamics at O(G^4)" 
by Z. Bern, J. Parra-Martinez, R. Roiban, M. S. Ruf, C.-H. Shen, M. P. Solon and M. Zeng.

M4cons = the complete G^4 contributions to the scattering.
This completes the calculation of https://arxiv.org/abs/2101.07254.

Ir4 = the G^4 term in the radial action for all conservative contributions, including the tail.

E[1]=Sqrt[m[1]^2+p^2]
E[2]=Sqrt[m[2]^2+p^2]
M = m[1]+m[2]
Etot = E[1]+E[2]
nu = m[1]m[2]/M^2
xi = E[1]E[2]/Etot^2
sigma = (E[1]E[2]+p^2)/(m[1]m[2])

H^{isotropic} = E_1 + E^2 + \sum_{n=1}^4 (G/r)^n * c_n
WARNING:  This Hamiltonian gives the correct scattering angle, 
          but beyond this it cannot be directly used for the bound
          orbits due to tail effect nonlocalities.

c1 = 1PM Hamiltonian coefficient
c2 = 2PM Hamiltonian coefficient
c3 = 3PM Hamiltonian coefficient
c4 = 4PM Hamiltonian coefficient

c4lower = the part of the 4PM Hamiltonian coefficient determined by
          c1, c2 and c3

M4p = test-body (Schwarzschild) part of the 3-loop amplitude
M4t = tail part of the 3-loop amplitude up to normalization
M4f = part of the first-order self-force regular in pInf=Sqrt[sigma^2-1]

The 2nd self force terms all cancel from the conservative contributions.

*)
 
c1 = (M^4*nu^2*(1 - 2*sigma^2))/(Etot^2*xi)

c2 = (M^5*nu^2*((3*(1 - 5*sigma^2))/4 - (M^3*nu^2*(1 - 2*sigma^2)^2*(1 - xi))/
        (2*Etot^3*xi^2) - (4*M*nu*sigma*(1 - 2*sigma^2))/(Etot*xi)))/(Etot^2*xi)

c3 = (M^6*nu^2*((-3*Etot*nu*(1 - 5*sigma^2)*(1 - 2*sigma^2))/
        (2*(1 + Etot/M)*M*(1 + sigma)) + (3 - 6*nu + 206*nu*sigma - 
         54*sigma^2 + 108*nu*sigma^2 + 4*nu*sigma^3)/12 + 
       (M^6*nu^4*(1 - 2*sigma^2)^3*(1 - 2*xi))/(2*Etot^6*xi^4) + 
       (2*M^4*nu^3*sigma*(1 - 2*sigma^2)^2*(3 - 4*xi))/(Etot^4*xi^3) - 
       (3*M*nu*sigma*(7 - 20*sigma^2))/(2*Etot*xi) - 
       (M^3*nu^2*(1 - 2*sigma^2)*(3 + (8*Etot)/M - 15*sigma^2 - 
          (80*Etot*sigma^2)/M - 3*xi + 15*sigma^2*xi))/(4*Etot^3*xi^2) - 
       (4*nu*(3 + 12*sigma^2 - 4*sigma^4)*ArcSinh[Sqrt[-1 + sigma]/Sqrt[2]])/
        Sqrt[-1 + sigma^2]))/(Etot^2*xi)

c4 = c4lower + (M^7*nu^2*(M4p + nu*(M4pi2+M4rem+4 M4t Log[Sqrt[sigma^2-1]/2])))/(4*Etot^2*xi)

c4lower = (M^8*nu^3*sigma*(-125 + 513*sigma^2))/(8*Etot^3*xi^2) + 
     (2*M^14*nu^7*sigma*(-1 + 2*sigma^2)^3*(5 - 12*xi + 3*xi^2))/
      (Etot^9*xi^6) - (M^16*nu^8*(1 - 2*sigma^2)^4*(5 - 15*xi + 6*xi^2 + 
        xi^3))/(8*Etot^11*xi^7) + (M^12*nu^6*(1 - 2*sigma^2)^2*
       (9*M*(-1 + 5*sigma^2)*(-1 + 2*xi) + 16*Etot*(-1 + 14*sigma^2)*
         (-2 + 3*xi)))/(8*Etot^8*xi^5) + 
     (M^8*nu^4*(-2880*Etot*M*sigma^2*(-1 + 2*sigma^2) + 
        288*Etot*M*(1 - 5*sigma^2)*(-1 + 6*sigma^2) + 
        216*Etot*M*(3 - 48*sigma^2 + 100*sigma^4) + 27*M^2*(1 - 5*sigma^2)^2*
         (-1 + xi) + (54*Etot^2*M^2*(1 - 5*sigma^2)^2*xi^2)/p^2 + 
        (48*Etot^2*M^2*(1 - 20*sigma^2 + 36*sigma^4)*xi^2)/p^2 + 
        (8*(3*M^3*(1 + sigma)^2*(-1 + sigma + 20*sigma^2 - 20*sigma^3 - 
             36*sigma^4 + 36*sigma^5)*(-1 + xi) + 8*Etot^3*(-37 - 103*sigma + 
             151*sigma^2 + 535*sigma^3 + 20*sigma^4 - 604*sigma^5 - 
             92*sigma^6 + 196*sigma^7)*xi + 3*Etot*M^2*(-1 + sigma)*
            (1 + sigma)^2*(-43 + sigma^2*(632 - 20*xi) + xi + 
             12*sigma^4*(-103 + 3*xi)) - 2*Etot^2*M*(1 + sigma)*
            (-63 + sigma^2*(981 - 643*xi) + 139*xi + 120*sigma*xi - 
             552*sigma^3*xi + 432*sigma^5*xi + 4*sigma^6*(450 + 29*xi) + 
             sigma^4*(-2718 + 256*xi))))/((Etot + M)*(-1 + sigma)*
          (1 + sigma)^2)))/(96*Etot^5*xi^3) - 
     (M^10*nu^5*(-1 + 2*sigma^2)*(4*Etot^4*M*(-15 + 200*sigma + 377*sigma^2 + 
          112*sigma^3 - 86*sigma^4)*xi^3 - 16*Etot^5*(3 - 25*sigma - 
          55*sigma^2 - 14*sigma^3 + 22*sigma^4)*xi^3 + 
        27*M^3*p^2*sigma*(-3 - 3*sigma + 10*sigma^2 + 10*sigma^3)*
         (-3 + 4*xi) + 2*Etot*M^2*p^2*(1 + sigma)*(-3*(-1 + xi)*xi + 
          54*sigma^2*(-1 + xi)*xi + 2*sigma^3*(-741 + 539*xi + xi^2) + 
          sigma*(387 - 427*xi + 103*xi^2)) + Etot^2*M*p^2*
         (6*(5 - 8*xi)*xi - 2*sigma^4*(1749 - 626*xi + 176*xi^2) + 
          sigma^3*(-3498 + 856*xi + 224*xi^2) + sigma*(819 - 724*xi + 
            400*xi^2) + sigma^2*(819 - 1078*xi + 880*xi^2)) - 
        4*Etot^3*(-(M^2*(-3 + 100*sigma + 157*sigma^2 + 56*sigma^3 + 
             2*sigma^4)*xi^3) + 2*p^2*(3*(-1 + xi)*xi + 
            sigma^2*(-36 + 55*xi - 55*xi^2) + sigma*(-36 + 25*xi - 25*xi^2) - 
            14*sigma^3*(-12 - xi + xi^2) + 2*sigma^4*(84 - 11*xi + 
              11*xi^2)))))/(12*Etot^6*(Etot + M)^2*p^2*(1 + sigma)*xi^4) + 
     ((4*M^8*nu^4*sigma*(15 - 118*sigma^2 + 132*sigma^4 - 40*sigma^6))/
        (Etot^3*(-1 + sigma^2)^(3/2)*xi^2) - 
       (4*M^10*nu^5*(3 + 6*sigma^2 - 28*sigma^4 + 8*sigma^6)*
         (p^2*(-1 + xi) + 2*Etot^2*xi^2))/(Etot^5*p^2*Sqrt[-1 + sigma^2]*
         xi^3))*ArcSinh[Sqrt[-1 + sigma]/Sqrt[2]]


(* G^4 contribution to radial action related to scattering angle via chi = - partial I_r/ partial J *)

Ir4 = - G^4 * M^7 * nu^2 * Pi * p^2/ (8 Etot J^3) * M4

(* Conservative scattering amplitude at order G^4*)
M4 = M4p + nu*(M4pi2+M4rem+4 M4t Log[Sqrt[sigma^2-1]/2])

(* Amplitude as normalized in equation (3) *)
M4cons = G^4 * M^7 * nu^2 * q * Pi^2 * M4  

(* Schwarzschild contribution *)
M4p = (-35*(1 - 18*sigma^2 + 33*sigma^4))/(8*(-1 + sigma^2))

(* tail *)
M4t=r[1]+r[2]*Log[(sigma+1)/2]+r[3]*ArcCosh[sigma]/Sqrt[sigma^2-1];

(* Piece that expands to terms proportional to Pi^2 in the near-static limit *)

M4pi2=r[4]*Pi^2+r[5]*EllipticK[(sigma-1)/(sigma+1)]EllipticE[(sigma-1)/(sigma+1)]+r[6]*EllipticK[(sigma-1)/(sigma+1)]^2+r[7]*EllipticE[(sigma-1)/(sigma+1)]^2;

(* Remainder *)

M4rem=r[8]+r[9]*Log[(sigma+1)/2]+r[10]*ArcCosh[sigma]/Sqrt[sigma^2-1]+r[11]*Log[sigma]+r[12]*Log[(sigma+1)/2]^2+r[13]*ArcCosh[sigma]/Sqrt[sigma^2-1]*Log[(sigma+1)/2]+
r[14]*ArcCosh[sigma]^2/(sigma^2-1)+
r[15]*PolyLog[2,(1-sigma)/2]+r[16]*PolyLog[2,(1-sigma)/(1+sigma) ]+r[17]*(PolyLog[2, -(Sqrt[-1 + sigma]/Sqrt[1 + sigma])] - 
    PolyLog[2, Sqrt[-1 + sigma]/Sqrt[1 + sigma]])/Sqrt[-1 + sigma^2] ;

(* reoccuring rational functions *)

g[1] = (2 + 15 sigma^2);
g[2] = sigma (19 + 15 sigma^2);
g[3] = (11 - 30 sigma^2 + 35 sigma^4);
g[4] = (-20 - 111 sigma^2 - 30 sigma^4 + 25 sigma^6);

(* rational function in tail *)

r[1] = 1/12*((1151 - 3336 sigma + 3148 sigma^2 - 912 sigma^3 + 
     339 sigma^4 - 552 sigma^5 + 210 sigma^6)/ (-1 + sigma^2));
r[2] = 4*g[1] - 2 g[2] - 1/2 g[3];
r[3] = (sigma (-3 + 2 sigma^2) )/(4 (-1 + sigma^2)) g[3];

(* rational function in finite part *)

r[4] = -1/3*g[4] - 16 /3 g[1] + 4/3 g[2];
r[5]= (-1183 - 2929 sigma - 2660 sigma^2 - 1200 sigma^3)/(
  2 (-1 + sigma) (1 + sigma));
r[6]= (834 + 2095 sigma + 1200 sigma^2)/(
  2 (-1 + sigma) (1 + sigma));
r[7]= (7 (169 + 380 sigma^2))/(4 (-1 + sigma));
r[8]= 1/(72 (-1 + sigma)^2 sigma^7 (1 + sigma)^2) (-45 + 207 sigma^2 - 
     1471 sigma^4 + 13349 sigma^6 - 38135 sigma^7 + 64424 sigma^8 - 
     32177 sigma^9 - 15056 sigma^10 - 25145 sigma^11 + 
     27952 sigma^12 + 33249 sigma^13 - 35360 sigma^14 + 
     4320 sigma^15 + 3600 sigma^16);
r[9]=-1/(3 (-1 + sigma^2))*
2 (-152 - 2520 sigma - 136 sigma^2 - 739 sigma^3 - 612 sigma^4 - 
       416 sigma^5 + 75 sigma^7);
r[10]= -2((sigma (-3 + 2 sigma^2))/( (-1 + sigma^2))) r[1] + 28*sigma^2 (-3 + 2 sigma^2);
r[11]=(
  4 sigma (-852 - 283 sigma^2 - 140 sigma^4 + 75 sigma^6))/(
  3 (-1 + sigma^2));
r[12]= 6 g[2] + g[3] - 1/2 g[4];
r[13]=-8 sigma (-3 + 2 sigma^2)/(-1 + sigma^2) g[1];
r[14]=-(((sigma^2) ((-3 + 2 sigma^2)^2)  )/(
    4 ((-1 + sigma^2)^2) )) g[3];
r[15]=-16 g[1] - 4 g[2] - g[4];
r[16]=-2 g[4];
r[17]=sigma (-3 + 2 sigma^2)/(-1 + sigma^2) (4 g[2] - 8 g[1]);

