Twist coefficients for the billiard map when the table is symmetric

0. Derivatives of radius along a curve

Clear[y,arc,radi,radi1,radi2,radi3,radi4];
y[x_,a_,b_,c_,d_]:=a*x^2+b*x^4+c*x^6+d*x^8; (*the parametric curve*)
arc[x_,a_,b_,c_,d_]:=Sqrt[1+D[y[x,a,b,c,d],x]^2];(*arc-length of the curve*)
radi[x_,a_,b_,c_,d_]:=arc[x,a,b,c,d]^3/(D[y[x,a,b,c,d],{x,2}]); (*radius of curvature*)
radi1[x_,a_,b_,c_,d_]:=D[radi[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D of arclength*)
radi2[x_,a_,b_,c_,d_]:=D[radi1[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D^2 of arclength*)
radi3[x_,a_,b_,c_,d_]:=D[radi2[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D^3 of arclength*)
radi4[x_,a_,b_,c_,d_]:=D[radi3[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D^4 of arclength*)

Simplify[{radi[x,a,b,c,d],radi1[x,a,b,c,d],radi2[x,a,b,c,d],radi3[x,a,b,c,d],radi4[x,a,b,c,d]}/.x->0,Element[{a,b,c,d},Reals]]

1. Derivatives of the billiard map (s,u)->((Subscript[s, 1]),(Subscript[u, 1]))

Clear[rad,tau,s1,u1];
(*set up the radius function*)
rad[0,r_,r2_,r4_]:=rad[0,r,r2,r4]=r;
Derivative[1,0,0,0][rad][0,r_,r2_,r4_]:=0;
Derivative[2,0,0,0][rad][0,r_,r2_,r4_]:=r2;
Derivative[3,0,0,0][rad][0,r_,r2_,r4_]:=0;
Derivative[4,0,0,0][rad][0,r_,r2_,r4_]:=r4;
(*set up the orbit length function and the billiard map*)tau[0,0,t_,r_,r2_,r4_]:=tau[0,0,t,r,r2,r4]=t;
u1[0,0,t_,r_,r2_,r4_]:=u1[0,0,t,r,r2,r4]=0;
s1[0,0,t_,r_,r2_,r4_]:=s1[0,0,t,r,r2,r4]=0;
Derivative[1,0,0,0,0,0][s1][s_,u_,t_,r_,r2_,r4_]:=tau[s,u,t,r,r2,r4]/(rad[s,r,r2,r4]*Sqrt[1-u1[s,u,t,r,r2,r4]^2])-Sqrt[1-u^2]/Sqrt[1-u1[s,u,t,r,r2,r4]^2];
Derivative[0,1,0,0,0,0][s1][s_,u_,t_,r_,r2_,r4_]:=(tau[s,u,t,r,r2,r4]/(Sqrt[1-u^2]*Sqrt[1-u1[s,u,t,r,r2,r4]^2]));
Derivative[1,0,0,0,0,0][u1][s_,u_,t_,r_,r2_,r4_]:=tau[s,u,t,r,r2,r4]/(rad[s,r,r2,r4]*rad[s1[s,u,t,r,r2,r4],r,r2,r4])-Sqrt[1-u^2]/rad[s1[s,u,t,r,r2,r4],r,r2,r4] -Sqrt[1-u1[s,u,t,r,r2,r4]^2]/rad[s,r,r2,r4];
Derivative[0,1,0,0,0,0][u1][s_,u_,t_,r_,r2_,r4_]:=tau[s,u,t,r,r2,r4]/(rad[s1[s,u,t,r,r2,r4],r,r2,r4]*Sqrt[1-u^2])-Sqrt[1-u1[s,u,t,r,r2,r4]^2]/Sqrt[1-u^2];
Derivative[1,0,0,0,0,0][tau][s_,u_,t_,r_,r2_,r4_]:=u-u1[s,u,t,r,r2,r4]*(tau[s,u,t,r,r2,r4]/(rad[s,r,r2,r4]*Sqrt[1-u1[s,u,t,r,r2,r4]^2])-Sqrt[1-u^2]/Sqrt[1-u1[s,u,t,r,r2,r4]^2]);
Derivative[0,1,0,0,0,0][tau][s_,u_,t_,r_,r2_,r4_]:=-u1[s,u,t,r,r2,r4]*(tau[s,u,t,r,r2,r4]/(Sqrt[1-u^2]*Sqrt[1-u1[s,u,t,r,r2,r4]^2]));

Higher-order derivatives of the billiard map (s,u)->((Subscript[s, 1]),(Subscript[u, 1]))

Clear[t,r,r2,r4,sderi,uderi];
t[v__]:=t[v]=v[[1]]; (*orbit length*)
r[v__]:=r[v]=v[[2]]; (*radius at the impact point*)
r2[v__]:=r2[v]=v[[3]]; (*second derivative of the radius at the impact point*)
r4[v__]:=r4[v]=v[[4]];(*fourth derivative of the radius at the impact point*)
sderi[j_,k_,v__]:=sderi[j,k,v]=(D[s1[svariable,uvariable,t[v],r[v],r2[v],r4[v]],{svariable,j},{uvariable,k}]/.svariable->0)/.uvariable->0;
uderi[j_,k_,v__]:=uderi[j,k,v]=(D[u1[svariable,uvariable,t[v],r[v],r2[v],r4[v]],{svariable,j},{uvariable,k}]/.svariable->0)/.uvariable->0;

2. Taylor expansion of the billiard map in the coordinate system (x,y)=(s/η, ηu)

The scaling and the eigenvalues

Clear[eta,lambda];
eta[v__]:=eta[v]=(-sderi[0,1,v]/uderi[1,0,v])^(1/4);
lambda[v__]:=lambda[v]=sderi[1,0,v]-I*Sqrt[-sderi[0,1,v]*uderi[1,0,v]];
(*since b10 = sinθ < 0 *)

Coefficients of (Subscript[x, 1])=∑(Subscript[a, jk])(x^j)(y^k), (Subscript[y, 1])=∑(Subscript[b, jk])(x^j)(y^k)

Clear[a10,a01,b10,b01,a30,a21,a12,a03,b30,b21,b12,b03,a50,a41,a32,a23,a14,a05,b50,b41,b32,b23,b14,b05];
a10[v__]:=a10[v]=sderi[1,0,v];a01[v__]:=a01[v]=eta[v]^(-2)*sderi[0,1,v];a30[v__]:=a30[v]=eta[v]^2*sderi[3,0,v]/6;a21[v__]:=a21[v]=sderi[2,1,v]/2;a12[v__]:=a12[v]=eta[v]^(-2)*sderi[1,2,v]/2;a03[v__]:=a03[v]=eta[v]^(-4)*sderi[0,3,v]/6;a50[v__]:=a50[v]=eta[v]^4*sderi[5,0,v]/5!;a41[v__]:=a41[v]=eta[v]^2*sderi[4,1,v]/24;a32[v__]:=a32[v]=sderi[3,2,v]/12;a23[v__]:=a23[v]=eta[v]^(-2)*sderi[2,3,v]/12;a14[v__]:=a14[v]=eta[v]^(-4)*sderi[1,4,v]/24;a05[v__]:=a05[v]=eta[v]^(-6)*sderi[0,5,v]/5!;
b10[v__]:=b10[v]=eta[v]^(2)*uderi[1,0,v];b01[v__]:=b01[v]=uderi[0,1,v];b30[v__]:=b30[v]=eta[v]^4*uderi[3,0,v]/6;b21[v__]:=b21[v]=eta[v]^2*uderi[2,1,v]/2;b12[v__]:=b12[v]=uderi[1,2,v]/2;b03[v__]:=b03[v]=eta[v]^(-2)*uderi[0,3,v]/6;b50[v__]:=b50[v]=eta[v]^6*uderi[5,0,v]/5!;b41[v__]:=b41[v]=eta[v]^4*uderi[4,1,v]/24;b32[v__]:=b32[v]=eta[v]^2*uderi[3,2,v]/12;b23[v__]:=b23[v]=uderi[2,3,v]/12;b14[v__]:=b14[v]=eta[v]^(-2)*uderi[1,4,v]/24;b05[v__]:=b05[v]=eta[v]^(-4)*uderi[0,5,v]/5!;

3. Coefficients of the coordinate transform (x,y)->(X,Y)=(∑(Subscript[p, jk])(x^j)(y^k),∑(Subscript[q, jk])(x^j)(y^k))

Clear[p30,p21,p12,p03,q30,q21,q12,q03,p50,p41,p32,p23,p14,p05,q50,q41,q32,q23,q14,q05];
p30[v__]:=p30[v]=-2 ((a30[v]+b21[v]-a12[v]-b03[v])a10[v]/8
+(b30[v]-a21[v]-b12[v]+a03[v])(a10[v]^2+a01[v] b10[v])/(16Sqrt[-a01[v] b10[v]]))+((a30[v]-b21[v]-a12[v]+b03[v])/(32a10[v])
-(b30[v]+a21[v]-b12[v]-a03[v])/(32Sqrt[-a01[v] b10[v]]));
p21[v__]:=p21[v]=3 ((a30[v]-b21[v]-a12[v]+b03[v])/(32Sqrt[-a01[v] b10[v]])+(b30[v]+a21[v]-b12[v]-a03[v])/(32a10[v]));
p12[v__]:=p12[v]=-6((a30[v]+b21[v]-a12[v]-b03[v])a10[v]/8
+(b30[v]-a21[v]-b12[v]+a03[v])(a10[v]^2+a01[v] b10[v])/(16Sqrt[-a01[v] b10[v]]))-3 ((a30[v]-b21[v]-a12[v]+b03[v])/(32a10[v])
-(b30[v]+a21[v]-b12[v]-a03[v])/(32Sqrt[-a01[v] b10[v]]));
p03[v__]:=p03[v]=4((b30[v]-a21[v]-b12[v]+a03[v])a10[v]/8-(a30[v]+b21[v]-a12[v]-b03[v])(a10[v]^2+a01[v] b10[v])/(16Sqrt[-a01[v] b10[v]])) -((a30[v]-b21[v]-a12[v]+b03[v])/(32Sqrt[-a01[v] b10[v]])+(b30[v]+a21[v]-b12[v]-a03[v])/(32a10[v]));
q30[v__]:=q30[v]=4 ((b30[v]-a21[v]-b12[v]+a03[v])a10[v]/8-(a30[v]+b21[v]-a12[v]-b03[v])(a10[v]^2+a01[v] b10[v])/(16Sqrt[-a01[v] b10[v]]))+((a30[v]-b21[v]-a12[v]+b03[v])/(32Sqrt[-a01[v] b10[v]])+(b30[v]+a21[v]-b12[v]-a03[v])/(32a10[v]));
q21[v__]:=q21[v]=6 ((a30[v]+b21[v]-a12[v]-b03[v])a10[v]/8
+(b30[v]-a21[v]-b12[v]+a03[v])(a10[v]^2+a01[v] b10[v])/(16Sqrt[-a01[v] b10[v]])) -3((a30[v]-b21[v]-a12[v]+b03[v])/(32a10[v])
-(b30[v]+a21[v]-b12[v]-a03[v])/(32Sqrt[-a01[v] b10[v]]));
q12[v__]:=q12[v]=-3 ((a30[v]-b21[v]-a12[v]+b03[v])/(32Sqrt[-a01[v] b10[v]])+(b30[v]+a21[v]-b12[v]-a03[v])/(32a10[v]));
q03[v__]:=q03[v]=2 ((a30[v]+b21[v]-a12[v]-b03[v])a10[v]/8
+(b30[v]-a21[v]-b12[v]+a03[v])(a10[v]^2+a01[v] b10[v])/(16Sqrt[-a01[v] b10[v]]))+((a30[v]-b21[v]-a12[v]+b03[v])/(32a10[v])
-(b30[v]+a21[v]-b12[v]-a03[v])/(32Sqrt[-a01[v] b10[v]]));
p50[v__]:=p50[v]=p21[v]*q30[v];
p41[v__]:=p41[v]=p21[v]*q21[v]+2p12[v]*q30[v];
p32[v__]:=p32[v]=p21[v]*q12[v]+2p12[v]*q21[v]+3p03[v]*q30[v];
p23[v__]:=p23[v]=p21[v]*q03[v]+2p12[v]*q12[v]+3p03[v]*q21[v];
p14[v__]:=p14[v]=2p12[v]*q03[v]+3p03[v]*q12[v];
p05[v__]:=p05[v]=3p03[v]*q03[v];q50[v__]:=q50[v]=q21[v]*q30[v];
q41[v__]:=q41[v]=q21[v]*q21[v]+2q12[v]*q30[v];
q32[v__]:=q32[v]=q21[v]*q12[v]+2q12[v]*q21[v]+3q03[v]*q30[v];
q23[v__]:=q23[v]=q21[v]*q03[v]+2q12[v]*q12[v]+3q03[v]*q21[v];
q14[v__]:=q14[v]=2q12[v]*q03[v]+3q03[v]*q12[v];
q05[v__]:=q05[v]=3q03[v]*q03[v];

4. Taylor expansion of the billiard map in terms of new coordinate system (X,Y): (X,Y)->((Subscript[X, 1]),(Subscript[Y, 1]))=(∑(Subscript[A, jk])(X^j)(Y^k),∑(Subscript[B, jk])(X^j)(Y^k))

Clear[aa30,aa21,aa12,aa03,bb30,bb21,bb12,bb03];
aa30[v__]:=aa30[v]=a30[v]+((b10[v])^3) p03[v]+a10[v] ((b10[v])^2) p12[v]+((a10[v])^2) b10[v]* p21[v]+((a10[v])^3) p30[v]-(a10[v]* p30[v]+a01[v]* q30[v]);
aa21[v__]:=aa21[v]=a21[v]+b01[v] (3 ((b10[v])^2) p03[v]+2 a10[v]* b10[v]* p12[v]+((a10[v])^2) p21[v])+a01[v] (((b10[v])^2) p12[v]+2 a10[v]* b10[v]* p21[v]+3 ((a10[v])^2) p30[v])-(a10[v]* p21[v]+a01[v]* q21[v]);aa12[v__]:=aa12[v]=a12[v]+((b01[v])^2) (3 b10[v]* p03[v]+a10[v]* p12[v])+2 a01[v]* b01[v] (b10[v]* p12[v]+a10[v]* p21[v])+((a01[v])^2) (b10[v]* p21[v]+3 a10[v]* p30[v])-(a10[v]* p12[v]+a01[v]* q12[v]);aa03[v__]:=aa03[v]=a03[v]+((b01[v])^3) p03[v]+a01[v] ((b01[v])^2) p12[v]+((a01[v])^2) b01[v]* p21[v]+((a01[v])^3) p30[v]-(a10[v]* p03[v]+a01[v]* q03[v]);
bb30[v__]:=bb30[v]=b30[v]+((b10[v])^3) q03[v]+a10[v] ((b10[v])^2) q12[v]+((a10[v])^2) b10[v]* q21[v]+((a10[v])^3) q30[v]-(b10[v]* p30[v]+b01[v]* q30[v]);
bb21[v__]:=bb21[v]=b21[v]+b01[v] (3 ((b10[v])^2) q03[v]+2 a10[v]* b10[v]* q12[v]+((a10[v])^2) q21[v])+a01[v] (((b10[v])^2) q12[v]+2 a10[v]* b10[v]* q21[v]+3 ((a10[v])^2) q30[v])-(b10[v]* p21[v]+b01[v]* q21[v]);bb12[v__]:=bb12[v]=b12[v]+((b01[v])^2) (3 b10[v]* q03[v]+a10[v]* q12[v])+2 a01[v]* b01[v] (b10[v]* q12[v]+a10[v]* q21[v])+((a01[v])^2) (b10[v]* q21[v]+3 a10[v]* q30[v])-(b10[v]* p12[v]+b01[v]* q12[v]);bb03[v__]:=bb03[v]=b03[v]+((b01[v])^3) q03[v]+a01[v] ((b01[v])^2) q12[v]+((a01[v])^2) b01[v]* q21[v]+((a01[v])^3) q30[v]-(b10[v]* p03[v]+b01[v]* q03[v]);

Clear[aa50,aa41,aa32,aa23,aa14,aa05];
aa50[v__]:=aa50[v]=a50[v]+3 ((b10[v])^2) b30[v] p03[v]+((b10[v])^5) p05[v]+(a30[v] ((b10[v])^2)+2 a10[v] b10[v] b30[v]) p12[v]+a10[v] ((b10[v])^4) p14[v]+(2 a10[v] a30[v] b10[v]+((a10[v])^2) b30[v]) p21[v]+((a10[v])^2) ((b10[v])^3) p23[v]+3 ((a10[v])^2) a30[v] p30[v]+((a10[v])^3) ((b10[v])^2) p32[v]+((a10[v])^4) b10[v] p41[v]+((a10[v])^5) p50[v]-(3 aa30[v] p30[v]+a10[v] p50[v]+aa21[v] q30[v]+a01[v] q50[v]);
aa41[v__]:=aa41[v]=a41[v]+((b10[v])^4) (5 b01[v] p05[v]+a01[v] p14[v])+2 a10[v] ((b10[v])^3) (2 b01[v] p14[v]+a01[v] p23[v])+((b10[v])^2) (3 b21[v] p03[v]+a21[v] p12[v]+3 ((a10[v])^2) (b01[v] p23[v]+a01[v] p32[v]))+2 b10[v] (a10[v] b21[v] p12[v]+a01[v] b30[v] p12[v]+a10[v] a21[v] p21[v]+a01[v] a30[v] p21[v]+b01[v] (3 b30[v] p03[v]+a30[v] p12[v]+((a10[v])^3) p32[v])+2 a01[v] ((a10[v])^3) p41[v])+a10[v] (a10[v] b21[v] p21[v]+2 a01[v] b30[v] p21[v]+3 a10[v] a21[v] p30[v]+6 a01[v] a30[v] p30[v]+b01[v] (2 b30[v] p12[v]+2 a30[v] p21[v]+((a10[v])^3) p41[v])+5 a01[v] ((a10[v])^3) p50[v])-(3 aa30[v] p21[v]+2 aa21[v] p30[v]+a10[v] p41[v]+aa21[v] q21[v]+2 aa12[v] q30[v]+a01[v] q41[v]);
aa32[v__]:=aa32[v]=a32[v]+3 ((b01[v])^2) b30[v] p03[v]+a30[v] ((b01[v])^2) p12[v]+2 a10[v] b01[v] b21[v] p12[v]+2 a01[v] b01[v] b30[v] p12[v]+2 a10[v] a21[v] b01[v] p21[v]+2 a01[v] a30[v] b01[v] p21[v]+((a10[v])^2) b12[v] p21[v]+2 a01[v] a10[v] b21[v] p21[v]+((a01[v])^2) b30[v] p21[v]+((b10[v])^3) (10 ((b01[v])^2) p05[v]+4 a01[v] b01[v] p14[v]+((a01[v])^2) p23[v])+3 ((a10[v])^2) a12[v] p30[v]+6 a01[v] a10[v] a21[v] p30[v]+3 ((a01[v])^2) a30[v] p30[v]+((a10[v])^3) ((b01[v])^2) p32[v]+((b10[v])^2) (3 b12[v] p03[v]+a12[v] p12[v]+6 a10[v] ((b01[v])^2) p14[v]+6 a01[v] a10[v] b01[v] p23[v]+3 ((a01[v])^2) a10[v] p32[v])+4 a01[v] ((a10[v])^3) b01[v] p41[v]+b10[v] (3 ((a10[v])^2) ((b01[v])^2) p23[v]+2 b01[v] (3 b21[v] p03[v]+a21[v] p12[v]+3 a01[v] ((a10[v])^2) p32[v])+2 (a10[v] b12[v] p12[v]+a01[v] b21[v] p12[v]+a10[v] a12[v] p21[v]+a01[v] a21[v] p21[v]+3 ((a01[v])^2) ((a10[v])^2) p41[v]))+10 ((a01[v])^2) ((a10[v])^3) p50[v]-(3 aa30[v] p12[v]+2 aa21[v] p21[v]+aa12[v] p30[v]+a10[v] p32[v]+aa21[v] q12[v]+2 aa12[v] q21[v]+3 aa03[v] q30[v]+a01[v] q32[v]);
aa23[v__]:=aa23[v]=a23[v]+6 b01[v] b10[v] b12[v] p03[v]+3 ((b01[v])^2) b21[v] p03[v]+10 ((b01[v])^3) ((b10[v])^2) p05[v]+a21[v] ((b01[v])^2) p12[v]+2 a12[v] b01[v] b10[v] p12[v]+a03[v] ((b10[v])^2) p12[v]+2 a10[v] b01[v] b12[v] p12[v]+2 a01[v] b10[v] b12[v] p12[v]+2 a01[v] b01[v] b21[v] p12[v]+4 a10[v] ((b01[v])^3) b10[v] p14[v]+6 a01[v] ((b01[v])^2) ((b10[v])^2) p14[v]+2 a10[v] a12[v] b01[v] p21[v]+2 a01[v] a21[v] b01[v] p21[v]+2 a03[v] a10[v] b10[v] p21[v]+2 a01[v] a12[v] b10[v] p21[v]+2 a01[v] a10[v] b12[v] p21[v]+((a01[v])^2) b21[v] p21[v]+b03[v] (3 ((b10[v])^2) p03[v]+2 a10[v] b10[v] p12[v]+((a10[v])^2) p21[v])+((a10[v])^2) ((b01[v])^3) p23[v]+6 a01[v] a10[v] ((b01[v])^2) b10[v] p23[v]+3 ((a01[v])^2) b01[v] ((b10[v])^2) p23[v]+3 a03[v] ((a10[v])^2) p30[v]+6 a01[v] a10[v] a12[v] p30[v]+3 ((a01[v])^2) a21[v] p30[v]+3 a01[v] ((a10[v])^2) ((b01[v])^2) p32[v]+6 ((a01[v])^2) a10[v] b01[v] b10[v] p32[v]+((a01[v])^3) ((b10[v])^2) p32[v]+6 ((a01[v])^2) ((a10[v])^2) b01[v] p41[v]+4 ((a01[v])^3) a10[v] b10[v] p41[v]+10 ((a01[v])^3) ((a10[v])^2) p50[v]-(3 aa30[v] p03[v]+2 aa21[v] p12[v]+aa12[v] p21[v]+a10[v] p23[v]+aa21[v] q03[v]+2 aa12[v] q12[v]+3 aa03[v] q21[v]+a01[v] q23[v]);
aa14[v__]:=aa14[v]=a14[v]+((b01[v])^4) (5 b10[v] p05[v]+a10[v] p14[v])+2 a01[v] ((b01[v])^3) (2 b10[v] p14[v]+a10[v] p23[v])+((b01[v])^2) (3 b12[v] p03[v]+a12[v] p12[v]+3 ((a01[v])^2) (b10[v] p23[v]+a10[v] p32[v]))+2 b01[v] (3 b03[v] b10[v] p03[v]+a10[v] b03[v] p12[v]+a03[v] b10[v] p12[v]+a01[v] b12[v] p12[v]+a03[v] a10[v] p21[v]+a01[v] a12[v] p21[v]+((a01[v])^3) b10[v] p32[v]+2 ((a01[v])^3) a10[v] p41[v])+a01[v] (2 a03[v] b10[v] p21[v]+a01[v] b12[v] p21[v]+2 b03[v] (b10[v] p12[v]+a10[v] p21[v])+6 a03[v] a10[v] p30[v]+3 a01[v] a12[v] p30[v]+((a01[v])^3) b10[v] p41[v]+5 ((a01[v])^3) a10[v] p50[v])-(2 aa21[v] p03[v]+aa12[v] p12[v]+a10[v] p14[v]+2 aa12[v] q03[v]+3 aa03[v] q12[v]+a01[v] q14[v]);
aa05[v__]:=aa05[v]=a05[v]+((b01[v])^5) p05[v]+a01[v] ((b01[v])^4) p14[v]+((a01[v])^2) ((b01[v])^3) p23[v]+((b01[v])^2) (3 b03[v] p03[v]+a03[v] p12[v]+((a01[v])^3) p32[v])+a01[v] b01[v] (2 b03[v] p12[v]+2 a03[v] p21[v]+((a01[v])^3) p41[v])+((a01[v])^2) (b03[v] p21[v]+3 a03[v] p30[v]+((a01[v])^3) p50[v])-(aa12[v] p03[v]+a10[v] p05[v]+3 aa03[v] q03[v]+a01[v] q05[v]);

Clear[bb50,bb41,bb32,bb23,bb14,bb05];
bb50[v__]:=bb50[v]=b50[v]+3 ((b10[v])^2) b30[v] q03[v]+((b10[v])^5) q05[v]+(a30[v] ((b10[v])^2)+2 a10[v] b10[v] b30[v]) q12[v]+a10[v] ((b10[v])^4) q14[v]+(2 a10[v] a30[v] b10[v]+((a10[v])^2) b30[v]) q21[v]+((a10[v])^2) ((b10[v])^3) q23[v]+3 ((a10[v])^2) a30[v] q30[v]+((a10[v])^3) ((b10[v])^2) q32[v]+((a10[v])^4) b10[v] q41[v]+((a10[v])^5) q50[v]-(3 bb30[v] p30[v]+b10[v] p50[v]+bb21[v] q30[v]+b01[v] q50[v]);
bb41[v__]:=bb41[v]=b41[v]+((b10[v])^4) (5 b01[v] q05[v]+a01[v] q14[v])+2 a10[v] ((b10[v])^3) (2 b01[v] q14[v]+a01[v] q23[v])+((b10[v])^2) (3 b21[v] q03[v]+a21[v] q12[v]+3 ((a10[v])^2) (b01[v] q23[v]+a01[v] q32[v]))+2 b10[v] (a10[v] b21[v] q12[v]+a01[v] b30[v] q12[v]+a10[v] a21[v] q21[v]+a01[v] a30[v] q21[v]+b01[v] (3 b30[v] q03[v]+a30[v] q12[v]+((a10[v])^3) q32[v])+2 a01[v] ((a10[v])^3) q41[v])+a10[v] (a10[v] b21[v] q21[v]+2 a01[v] b30[v] q21[v]+3 a10[v] a21[v] q30[v]+6 a01[v] a30[v] q30[v]+b01[v] (2 b30[v] q12[v]+2 a30[v] q21[v]+((a10[v])^3) q41[v])+5 a01[v] ((a10[v])^3) q50[v])-(3 bb30[v] p21[v]+2 bb21[v] p30[v]+b10[v] p41[v]+bb21[v] q21[v]+2 bb12[v] q30[v]+b01[v] q41[v]);
bb32[v__]:=bb32[v]=b32[v]+3 ((b01[v])^2) b30[v] q03[v]+a30[v] ((b01[v])^2) q12[v]+2 a10[v] b01[v] b21[v] q12[v]+2 a01[v] b01[v] b30[v] q12[v]+2 a10[v] a21[v] b01[v] q21[v]+2 a01[v] a30[v] b01[v] q21[v]+((a10[v])^2) b12[v] q21[v]+2 a01[v] a10[v] b21[v] q21[v]+((a01[v])^2) b30[v] q21[v]+((b10[v])^3) (10 ((b01[v])^2) q05[v]+4 a01[v] b01[v] q14[v]+((a01[v])^2) q23[v])+3 ((a10[v])^2) a12[v] q30[v]+6 a01[v] a10[v] a21[v] q30[v]+3 ((a01[v])^2) a30[v] q30[v]+((a10[v])^3) ((b01[v])^2) q32[v]+((b10[v])^2) (3 b12[v] q03[v]+a12[v] q12[v]+6 a10[v] ((b01[v])^2) q14[v]+6 a01[v] a10[v] b01[v] q23[v]+3 ((a01[v])^2) a10[v] q32[v])+4 a01[v] ((a10[v])^3) b01[v] q41[v]+b10[v] (3 ((a10[v])^2) ((b01[v])^2) q23[v]+2 b01[v] (3 b21[v] q03[v]+a21[v] q12[v]+3 a01[v] ((a10[v])^2) q32[v])+2 (a10[v] b12[v] q12[v]+a01[v] b21[v] q12[v]+a10[v] a12[v] q21[v]+a01[v] a21[v] q21[v]+3 ((a01[v])^2) ((a10[v])^2) q41[v]))+10 ((a01[v])^2) ((a10[v])^3) q50[v]-(3 bb30[v] p12[v]+2 bb21[v] p21[v]+bb12[v] p30[v]+b10[v] p32[v]+bb21[v] q12[v]+2 bb12[v] q21[v]+3 bb03[v] q30[v]+b01[v] q32[v]);
bb23[v__]:=bb23[v]=b23[v]+6 b01[v] b10[v] b12[v] q03[v]+3 ((b01[v])^2) b21[v] q03[v]+10 ((b01[v])^3) ((b10[v])^2) q05[v]+a21[v] ((b01[v])^2) q12[v]+2 a12[v] b01[v] b10[v] q12[v]+a03[v] ((b10[v])^2) q12[v]+2 a10[v] b01[v] b12[v] q12[v]+2 a01[v] b10[v] b12[v] q12[v]+2 a01[v] b01[v] b21[v] q12[v]+4 a10[v] ((b01[v])^3) b10[v] q14[v]+6 a01[v] ((b01[v])^2) ((b10[v])^2) q14[v]+2 a10[v] a12[v] b01[v] q21[v]+2 a01[v] a21[v] b01[v] q21[v]+2 a03[v] a10[v] b10[v] q21[v]+2 a01[v] a12[v] b10[v] q21[v]+2 a01[v] a10[v] b12[v] q21[v]+((a01[v])^2) b21[v] q21[v]+b03[v] (3 ((b10[v])^2) q03[v]+2 a10[v] b10[v] q12[v]+((a10[v])^2) q21[v])+((a10[v])^2) ((b01[v])^3) q23[v]+6 a01[v] a10[v] ((b01[v])^2) b10[v] q23[v]+3 ((a01[v])^2) b01[v] ((b10[v])^2) q23[v]+3 a03[v] ((a10[v])^2) q30[v]+6 a01[v] a10[v] a12[v] q30[v]+3 ((a01[v])^2) a21[v] q30[v]+3 a01[v] ((a10[v])^2) ((b01[v])^2) q32[v]+6 ((a01[v])^2) a10[v] b01[v] b10[v] q32[v]+((a01[v])^3) ((b10[v])^2) q32[v]+6 ((a01[v])^2) ((a10[v])^2) b01[v] q41[v]+4 ((a01[v])^3) a10[v] b10[v] q41[v]+10 ((a01[v])^3) ((a10[v])^2) q50[v]-(3 bb30[v] p03[v]+2 bb21[v] p12[v]+bb12[v] p21[v]+b10[v] p23[v]+bb21[v] q03[v]+2 bb12[v] q12[v]+3 bb03[v] q21[v]+b01[v] q23[v]);
bb14[v__]:=bb14[v]=b14[v]+((b01[v])^4) (5 b10[v] q05[v]+a10[v] q14[v])+2 a01[v] ((b01[v])^3) (2 b10[v] q14[v]+a10[v] q23[v])+((b01[v])^2) (3 b12[v] q03[v]+a12[v] q12[v]+3 ((a01[v])^2) (b10[v] q23[v]+a10[v] q32[v]))+2 b01[v] (3 b03[v] b10[v] q03[v]+a10[v] b03[v] q12[v]+a03[v] b10[v] q12[v]+a01[v] b12[v] q12[v]+a03[v] a10[v] q21[v]+a01[v] a12[v] q21[v]+((a01[v])^3) b10[v] q32[v]+2 ((a01[v])^3) a10[v] q41[v])+a01[v] (2 a03[v] b10[v] q21[v]+a01[v] b12[v] q21[v]+2 b03[v] (b10[v] q12[v]+a10[v] q21[v])+6 a03[v] a10[v] q30[v]+3 a01[v] a12[v] q30[v]+((a01[v])^3) b10[v] q41[v]+5 ((a01[v])^3) a10[v] q50[v])-(2 bb21[v] p03[v]+bb12[v] p12[v]+b10[v] p14[v]+2 bb12[v] q03[v]+3 bb03[v] q12[v]+b01[v] q14[v]);
bb05[v__]:=bb05[v]=b05[v]+3 ((b01[v])^2) b03[v] q03[v]+((b01[v])^5) q05[v]+(a03[v] ((b01[v])^2)+2 a01[v] b01[v] b03[v]) q12[v]+a01[v] ((b01[v])^4) q14[v]+(2 a01[v] a03[v] b01[v]+((a01[v])^2) b03[v]) q21[v]+((a01[v])^2) ((b01[v])^3) q23[v]+3 ((a01[v])^2) a03[v] q30[v]+((a01[v])^3) ((b01[v])^2) q32[v]+((a01[v])^4) b01[v] q41[v]+((a01[v])^5) q50[v]-(bb12[v] p03[v]+b10[v] p05[v]+3 bb03[v] q03[v]+b01[v] q05[v]);

5. Twist Coefficients (Subscript[τ, 1]) and (Subscript[τ, 2])

Clear[twist1,twist2];
twist1[v__]:=twist1[v]=(a10[v](3b30[v]-a21[v]+b12[v]-3a03[v])+Sqrt[-a01[v]*b10[v]](3a30[v]+b21[v]+a12[v]+3b03[v]))/8;
twist2[v__]:=twist2[v]=(a10[v](-2aa41[v]-2aa23[v]-10aa05[v]+10bb50[v]+2bb32[v]+2bb14[v])+Sqrt[-a01[v]*b10[v]](10aa50[v]+2aa32[v]+2aa14[v]+2bb41[v]+2bb23[v]+10bb05[v]))/32;

Simplified formulas

Clear[tau1, tau2];
tau1[L_,R_,R2_]:=(1/(8R))-(L/(8(2R-L)))R2;
tau2[L_,R_,R2_,R4_]:=(1/64)(((3(7(R^2)-8R L+2(L^2)))/(4(R^2)(R-L)(((2R-L))^(1/2))(L^(1/2))))-(((L^(1/2))(27(R^2)-40R L+10(L^2)))/(6R(R-L)(((2R-L))^(3/2))))R2+(((L^(3/2))(31(R^2)-36R L+6(L^2)))/(12(R-L)(((2R-L))^(5/2))))(((R2))^2)-(((L^(3/2))R)/(3(((2R-L))^(3/2))))R4);

Integrability test: locally analytically integrable iff convergent Birkhoff Normal form. At (λ^4)=1: it is necessary to have (Subscript[c, 03])=0

Clear[laiTest];
laiTest[L_,R_,R2_]:=Simplify[a30[{L,R,R2}]+I b30[{L,R,R2}]+I a21[{L,R,R2}]−b21[{L,R,R2}]−a12[{L,R,R2}]−I b12[{L,R,R2}]−I a03[{L,R,R2}]+b03[{L,R,R2}],R>L&&L>0]