Birkhoff normal form and twist coefficients of periodic orbits of billiards​
By Pengfei Zhang and Xin Jin

Article abstract

In this paper we study the Birkhoff Normal Form around elliptic periodic points for a variety of dynamical billiards. We give an explicit construction of the Birkhoff transformation and obtain explicit formulas for the first two twist coefficients in terms of the geometric parameters of the billiard table. As an application, we obtain characterizations of the nonlinear stability and local analytic integrability of the billiards around the elliptic periodic points.

Introduction

Part 1 Symmetric billiards can be viewed as a special case of Part 2 Asymmetric billiards. The functions in Part 1 are redefined in Part 2. Evaluating the whole notebook would allow functions in Part 2 to override the functions in part 1.
To find the twist coefficients of symmetric billiards, please only evaluate Part 1. This part is used in Section 3 and Section 4 of the aforementioned paper.
To find the twist coefficients of asymmetric billiards, please only evaluate Part 2. This part is used in Section 5 of the aforementioned paper.
Part 3 is independent Part 1 and Part 2. It shows how to find the unknown polynomials under the Ansatz of formulation of
τ
2
. See Equation (5.12) Section 5 of the aforementioned paper.
Part 1. Twist coefficients for symmetric billiards

0. Derivatives of radius along a curve

In[]:=
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;(*theparametriccurve*)​​arc[x_,a_,b_,c_,d_]:=Sqrt[1+D[y[x,a,b,c,d],x]^2];(*arc-lengthofthecurve*)​​radi[x_,a_,b_,c_,d_]:=arc[x,a,b,c,d]^3/(D[y[x,a,b,c,d],{x,2}]);(*radiusofcurvature*)​​radi1[x_,a_,b_,c_,d_]:=D[radi[x,a,b,c,d],x]/arc[x,a,b,c,d];(*Dofarclength*)​​radi2[x_,a_,b_,c_,d_]:=D[radi1[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D^2ofarclength*)​​radi3[x_,a_,b_,c_,d_]:=D[radi2[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D^3ofarclength*)​​radi4[x_,a_,b_,c_,d_]:=D[radi3[x,a,b,c,d],x]/arc[x,a,b,c,d];(*D^4ofarclength*)
In[]:=
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]]
Out[]=

1
2a
,0,6a-
6b
2
a
,0,-
12(2
6
a
+4
3
a
b-36
2
b
+15ac)
3
a


1. Derivatives of the billiard map
(s,u)(
s
1
,
u
1
)

In[]:=
Clear[rad,tau,s1,u1];​​(*setuptheradiusfunction*)​​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;​​(*setuptheorbitlengthfunctionandthebilliardmap*)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)(
s
1
,
u
1
)
In[]:=
Clear[t,r,r2,r4,sderi,uderi];​​t[v__]:=t[v]=v[[1]];(*orbitlength*)​​r[v__]:=r[v]=v[[2]];(*radiusattheimpactpoint*)​​r2[v__]:=r2[v]=v[[3]];(*secondderivativeoftheradiusattheimpactpoint*)​​r4[v__]:=r4[v]=v[[4]];(*fourthderivativeoftheradiusattheimpactpoint*)​​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
In[]:=
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]];​​(*sinceb10=sinθ<0*)
Coefficients of
x
1
=
a
jk
j
x
k
y
,
y
1
=
b
jk
j
x
k
y
In[]:=
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)=
p
jk
j
x
k
y
,
q
jk
j
x
k
y

In[]:=

4. Taylor expansion of the billiard map in terms of new coordinate system
(X,Y)
:


5. Twist coefficients
τ
1
and
τ
2

To find
τ
1
and
τ
2
, use the input of the form
v={L,R,
R
2
,
R
4
}
, where L is the the length, R is the radius,
​
R
2
is the second - derivative, and
R
4
is the four - derivative
In[]:=
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
In[]:=
Clear[tau1,tau2];​​tau1[L_,R_,R2_]:=
1
8R
-
L
8(2R-L)
R2;​​tau2[L_,R_,R2_,R4_]:=
1
64
3(7
2
R
-8RL+2
2
L
)
4
2
R
(R-L)
1/2
(2R-L)
1/2
L
-
1/2
L
(27
2
R
-40RL+10
2
L
)
6R(R-L)
3/2
(2R-L)
R2+
3/2
L
(31
2
R
-36RL+6
2
L
)
12(R-L)
5/2
(2R-L)
2
(R2)
-
3/2
L
R
3
3/2
(2R-L)
R4;
Integrability test: locally analytically integrable iff convergent Birkhoff Normal form. At
4
λ
=1
: it is necessary to have
c
03
=0
In[]:=
Clear[laiTest];​​laiTest[L_,R_,R2_]:=Simplify[a30[{L,R,R2}]+Ib30[{L,R,R2}]+Ia21[{L,R,R2}]−b21[{L,R,R2}]−a12[{L,R,R2}]−Ib12[{L,R,R2}]−Ia03[{L,R,R2}]+b03[{L,R,R2}],R>L&&L>0]
Part 2. Twist coefficients for asymmetric billiards

1. Derivatives of the billiard map

(s,u)(
s
1
,
u
1
)(
s
2
,
u
2
)
Functions for the map from point A to point B
Clear[rad0,rad1,tau,s1,u1];​​(*setuptheradiusfunction*)​​rad0[0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rad0[0,t,ra,ra2,ra4,rb,rb2,rb4]=ra;​​Derivative[1,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;​​Derivative[2,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=ra2;​​Derivative[3,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;​​Derivative[4,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=ra4;​​rad1[0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rad1[0,t,ra,ra2,ra4,rb,rb2,rb4]=rb;​​Derivative[1,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;​​Derivative[2,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rb2;​​Derivative[3,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;​​Derivative[4,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rb4;​​(*setuptheorbitlengthfunctionandthebilliardmap*)tau[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=t;​​u1[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=u1[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=0;​​s1[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=s1[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=0;​​Derivative[1,0,0,0,0,0,0,0,0][s1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])-
Sqrt[1-u^2]
Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]
;​​Derivative[0,1,0,0,0,0,0,0,0][s1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=
tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]
Sqrt[1-u^2]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]
;​​Derivative[1,0,0,0,0,0,0,0,0][u1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]*rad1[s1[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4])-Sqrt[1-u^2]/rad1[s1[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4]-
Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]
rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]
;​​Derivative[0,1,0,0,0,0,0,0,0][u1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad1[s1[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u^2])-
Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]
Sqrt[1-u^2]
;​​Derivative[1,0,0,0,0,0,0,0,0][tau][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=u-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]*tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4](rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])-
Sqrt[1-u^2]
Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]
;​​Derivative[0,1,0,0,0,0,0,0,0][tau][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]*
tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]
Sqrt[1-u^2]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]
;
Functions for the map from point B to point A
Evaluate along the orbit from A to B and then back to A

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

3. Coefficients of the coordinate transform

5. Twist Coefficients

3.1 Find the symmetric part symP of P

4.1. The symmetric part of S(L,R,S)

4.2. The antisymmetric part of S