# Data-driven Discovery of The Quadrotor Equations of Motion Via Sparse
Identification of Nonlinear Dynamics

Zeyad M. Manaa
Aerospace Engineering Department
King Fahd University of Petroleum and Minerals
Dhahran, 31261, Saudi Arabian
&Mohammed R. Elbalshy
Aerospace Engineering Department
King Fahd University of Petroleum and Minerals
Dhahran, 31261, Saudi Arabian
&Ayman M. Abdallah
Aerospace Engineering Department
King Fahd University of Petroleum and Minerals
Dhahran, 31261, Saudi Arabian
Corresponding author: Z. M. M<EMAIL_ADDRESS>
Author contributions: Z. M. M formulated research; Z. M. M and M. R. E
performed research and simulations; Z. M. M, M. R. E, and A. M. A analyzed
data; M. R. E and Z. M. M wrote the manuscript; and Z. M. M, M. R. E, and A.
M. A reviewed and finalized the manuscript.
Z. M. M would like to acknowledge the support provided by the Deanship of
Research Oversight and Coordination at King Fahd University of Petroleum and
Minerals (KFUPM) under Research Grant xxxxx.

###### Abstract

Dynamical systems provide a mathematical framework for understanding complex
physical phenomena. The mathematical formulation of these systems plays a
crucial role in numerous applications; however, it often proves to be quite
intricate. Fortunately, data can be readily available through sensor
measurements or numerical simulations. In this study, we employ the Sparse
Identification of Nonlinear Dynamics (SINDy) algorithm to extract a
mathematical model solely from data. The influence of the hyperparameter
$\lambda$ on the sparsity of the identified dynamics is discussed.
Additionally, we investigate the impact of data size and the time step between
snapshots on the discovered model. To serve as a data source, a ground truth
mathematical model was derived from the first principals, we focus on modeling
the dynamics of a generic 6 Degrees of Freedom (DOF) quadrotor. For the scope
of this initial manuscript and for simplicity and algorithm validation
purposes, we specifically consider a sub-case of the 6 DOF system for
simulation, restricting the quadrotor’s motion to a 2-dimensional plane (i.e.
3 DOF). To evaluate the efficacy of the SINDy algorithm, we simulate three
cases employing a Proportional-Derivative (PD) controller for the 3 DOF case
including different trajectories. The performance of SINDy model is assessed
through the evaluation of absolute error metrics and root mean squared error
(RMSE). Interestingly, the predicted states exhibit at most a RMSE of order of
magnitude approximately $10^{-4}$, manifestation of the algorithm’s
effectiveness. This research highlights the application of the SINDy algorithm
in extracting the quadrotor mathematical models from data. We also try to
investigate the effect of noisy measurements on the algorithm efficacy. The
successful modeling of the 3 DOF quadrotor dynamics demonstrates the
algorithm’s potential, while the evaluation metrics validate its performance,
thereby clearing the way for more applications in the realm of unmanned aerial
vehicles.

_K_ eywords Keywords: dynamical systems $\cdot$ machine learning $\cdot$
sparse regression $\cdot$ system identification $\cdot$ quadrotor $\cdot$
numerical simulations $\cdot$ optimization

##  1 Introduction

Dynamical systems provide a mathematical representation for describing the
world. Formally, dynamical systems are concerned with the analysis, and
interpretation of the behavior of sets of differential equations that trace
the advancement of a system’s state across time [1, 2]. Classical dynamics has
been discussed for years and applied to a wide range of applications; it is
far too complex to be stated in a few words. However, the combination of big
data, machine learning techniques, and statistical learning is driving a
revolution in data-driven dynamical modeling and control of complex systems,
with analytical derivations being replaced by data-driven methodology [3]. In
most real-world circumstances, data is so abundant that it cannot be
comprehended. Furthermore, the physical principles that control these data are
frequently complex; this is true for most physical concerns, such as climate
science, epidemiology, and finance, to name a few. Researchers [4, 5, 6]
attempt to decompose data in order to acquire insights into these massive, yet
unexplained, datasets. As a result, the abundance of data remains a challenge
because it is inherently imperfect. Consequently, the rise of data-driven
dynamics is paving the way for such difficulties.

The outcome of this development in data science was not initially noticeable
on dynamical systems [7], but recently a lot of research has been geared
towards that area. Bongard and Lipson [8] and Schmidt and Lipson [9] was able
to glean information about a nonlinear dynamical system’s structure which was
used after that to find the nonlinear differential equations [10] by symbolic
regression. Dynamic Mode Decomposition (DMD) was first introduced to fluid
dynamics by Schmid [11]. He extracted spatiotemporal structures from high
dimensional data using DMD to analyse it. Researchers have extended the DMD
method to explore other interesting intersections, including but not limited
to the Koopman operator [12], extended DMD [13] to allow for approximation for
Koopman operator, kernel DMD [14, 15], time-delay DMD [16, 17]. Moreover,
Proctor et al. [18] extended the DMD to account for control input. These
methods allow for linear system identification. Although these techniques rely
on a precise set of coordinates for dynamics linearization, they are
beneficial for developing dynamics-based linear models that progress high-
dimensional observations across time [19, 20] which are useful for control
purposes [21]. Recent studies [22, 23] have examined the use of deep learning
techniques to choose the proper set of coordinates to be used in both DMD and
extended DMD.

While linear dynamical systems have several great qualities, there are many
complex and intriguing dynamical phenomena that a linear model cannot properly
capture. This encourages ways for discovering models of nonlinear dynamical
systems. Consequently, another breakthrough on data-driven dynamics happened
after Brunton et al. [7] revealed the method of discovering nonlinear
equations of motion using SINDy algorithm. They combined compressed sensing
[24, 25, 26] with sparse regression [27, 28] and came up with the new
technique of SINDy, which is based on the idea that most dynamical systems
contain a small number of active terms that account for the majority of
dynamics. Recently, SINDy algorithm has been generalized to include control
inputs [3], to include tensor bases [29], to discover partial differential
equation [30, 31], and to account for noisy data as a kind of challenge [32,
33]. The usage of SINDy algorithm thereafter begin to be a workhorse in the
fields of optics, chemical engineering, robotics, and disease modeling [34,
35, 36, 37, 38].

Despite the increasing body of literature in the field, there has been a
dearth of research focused on the domain of aerospace, especially the Unmanned
Aerial Vehicles (UAV). A few studies discussed data-driven dynamics and
control as an application of the UAVs. Manzoor et al. [39] proposed a novel
approach for controlling ducted fan aerial vehicles (DFAVs) by integrating
model predictive control (MPC) and physics-informed machine learning. They
combined the physics-based model with the data-driven model to account for
sudden changes in the system dynamics motivated by [40, 41]. Also, Kaiser et
al. [3] proposed an integration between SINDy and MPC. They found that the
SINDy method needs less data compared to other popular machine learning
approaches, such as neural networks. This reduces the likelihood of over-
fitting. In contrast to [39], when combined with MPC, the SINDy algorithm
yields models that are both computationally feasible and precise, even when
trained on limited data.

Motivated by this, SINDy algorithm will be utilized to discover the EOM of
quadrotor using numerical simulation data at first, and for the next phase,
experimental data might be used. Moreover, a study on SINDy’s sparsity
hyperparameter effect on model discovery will be examined along with the
effect of the sampling time of data collected whether numerically or
experimentally.

##  2 Modelling and Simulation

###  2.1 Coordinate Systems

The coordinate systems for the quadrotor are shown in Fig. 1. The inertial
reference frame, $\mathcal{I}$, is defined by

$\mathcal{I}=\\{\mathcal{O},\boldsymbol{x}_{\mathcal{I}},\boldsymbol{y}_{\mathcal{I}},\boldsymbol{z}_{\mathcal{I}}\\}$
(1)

where, $\mathcal{O}$ is the frame’s origin, $\mathbf{x}_{\mathcal{I}}$,
$\mathbf{y}_{\mathcal{I}}$, $\mathbf{z}_{\mathcal{I}}$ are unit vectors
defining the frame axes. The body frame, $\mathcal{B}$, is fixed to the center
of mass of the quadrotor defined by

$\mathcal{B}=\\{\mathcal{C},\boldsymbol{x}_{\mathcal{B}},\boldsymbol{y}_{\mathcal{B}},\boldsymbol{z}_{\mathcal{B}}\\}$
(2)

where, $\mathcal{C}$ is the frame’s origin, $\mathbf{x}_{\mathcal{B}}$,
$\mathbf{y}_{\mathcal{B}}$, $\mathbf{z}_{\mathcal{B}}$ are unit vectors
defining the frame axes.

Figure 1: Quadrotor with the body and the inertial reference frames.

We may acquire the transformation matrix that transition between the body
frame to the inertial frame by using the $\psi-\theta-\phi$ sequence, which
indicates the yaw, pitch and roll angles respectively.

${}^{\mathcal{I}}\mathcal{R}_{\mathcal{B}}=\left[\begin{array}[]{ccc}c\psi
c\theta-s\phi s\psi s\theta&-c\phi s\psi&c\psi s\theta+c\theta s\phi s\psi\\\
c\theta s\psi+c\psi s\phi s\theta&c\phi c\psi&s\psi s\theta-c\psi c\theta
s\phi\\\ -c\phi s\theta&s\phi&c\phi c\theta\end{array}\right]$ (3)

Giving the fact that this matrix is orthonormal

${}^{\mathcal{I}}\mathcal{R}_{\mathcal{B}}{}^{\mathcal{I}}\mathcal{R}^{-1}_{\mathcal{B}}={}^{\mathcal{I}}\mathcal{R}_{\mathcal{B}}{}^{\mathcal{I}}\mathcal{R}^{\top}_{\mathcal{B}}=\texttt{eye(3)}$
(4)

where eye(3) is a $3\times 3$ identity matrix.

###  2.2 Quadrotor Dynamics

We assume the quadrotor is a rigid body with 6 DOF, having a mass m and
Intertia matrix $\boldsymbol{J}=\texttt{diag}(J_{x},J_{y},J_{z})$. Let
$\boldsymbol{r}$ denotes the position of $\mathcal{C}$ in $\mathcal{I}$, the
$\sum\boldsymbol{F}$ denotes the summation of forces acting upon the body and
$T_{n}$ denotes the motor force. The equations governing the motion of
$\mathcal{C}$ is derived by applying Newton’s law to the translational motion
considering the control actuation applied on the body coordinate.

$\displaystyle m\ddot{\boldsymbol{r}}^{\mathcal{I}}$
$\displaystyle=m\left(\ddot{\boldsymbol{r}}^{\mathcal{B}}+{\boldsymbol{\omega}^{\mathcal{B/I}}\times{\dot{\boldsymbol{r}}}^{\mathcal{B}}}\right)=\sum{\boldsymbol{F}}$
(5) $\displaystyle\sum{\boldsymbol{F}}$
$\displaystyle={\boldsymbol{F}_{g}}+{\boldsymbol{F}}_{\text{th}}=\begin{bmatrix}0\\\
0\\\
-mg\end{bmatrix}+{}^{\mathcal{I}}\mathcal{R}_{\mathcal{B}}\begin{bmatrix}0\\\
0\\\ \sum_{n=1}^{4}{T_{n}}\end{bmatrix}$
$\displaystyle\ddot{\boldsymbol{r}}^{\mathcal{B}}$
$\displaystyle=\frac{1}{m}\sum{\boldsymbol{F}}-{\boldsymbol{\omega}^{\mathcal{B/I}}\times{\dot{\boldsymbol{r}}}^{\mathcal{B}}}$

We also employ Euler’s equation to model the attitude dynamics such that $J$
is the quadrotor’s moment of inertia, $\boldsymbol{\omega}$ is the angular
velocity, $\boldsymbol{h}$ is the angular momentum defined as
$\boldsymbol{h}=\boldsymbol{J}\boldsymbol{\omega}$, the $\sum\boldsymbol{M}$
is the summation of the moments acting upon the system.

$\dot{\boldsymbol{\omega}}=\boldsymbol{J}^{-1}\left({\sum\boldsymbol{M}}-\boldsymbol{\omega}\times{\boldsymbol{J}\boldsymbol{\omega}}\right)$
(6)

The relation between quadrotor velocity in the body frame and position in the
inertial frame can be expressed as
$\dot{\boldsymbol{r}}^{\mathcal{I}}={}^{\mathcal{I}}\mathcal{R}_{\mathcal{B}}\dot{\boldsymbol{r}}^{\mathcal{B}}$.
Also, to get the attitude of the quadrotor we define the relation between
quadrotor compenets of angular velocity in the body frame and the roll, pitch
and yaw derivatives in the inertial frame by
$\boldsymbol{\omega}=\mathcal{T}\boldsymbol{\dot{a}}$. In summary, the overall
dynamics of the system in hand are given by equation 7.

$\displaystyle\dot{\boldsymbol{r}}^{\mathcal{I}}$
$\displaystyle={}^{\mathcal{I}}\mathcal{R}_{\mathcal{B}}\dot{\boldsymbol{r}}^{\mathcal{B}}$
(7) $\displaystyle\ddot{\boldsymbol{r}}^{\mathcal{B}}$
$\displaystyle=\frac{1}{m}\sum{\boldsymbol{F}}-{\boldsymbol{\omega}^{\mathcal{B/I}}\times{\dot{\boldsymbol{r}}}^{\mathcal{B}}}$
$\displaystyle\boldsymbol{\omega}$
$\displaystyle=\mathcal{T}\boldsymbol{\dot{a}}$
$\displaystyle\dot{\boldsymbol{\omega}}$
$\displaystyle=\boldsymbol{J}^{-1}\left({\boldsymbol{M}}-\boldsymbol{\omega}\times{\boldsymbol{J}\boldsymbol{\omega}}\right)$

###  2.3 Control inputs

The first control actuation input is defined as the sum of the forces from
each of the four rotors. So,

$u_{1}=[T_{1}+T_{2}+T_{3}+T_{4}]$ (8)

As shown in Fig 1, rotors one and three both rotate in the negative
${\boldsymbol{z}_{\mathcal{B}}}$ direction. Rotors two and four both rotate in
the positive ${\boldsymbol{z}_{\mathcal{B}}}$ direction. Since the moment
exerted on the quadrotor oppose the direction of blade rotation. So, the
moments ${M_{1}}$ and ${M_{2}}$ are directed to the positive
${\boldsymbol{z}_{\mathcal{B}}}$ direction, while the moments ${M_{3}}$ and
${M_{4}}$ are directed to the positive ${\boldsymbol{z}_{\mathcal{B}}}$
direction. Recalling to Eq. 6, it can be rewritten as

${\boldsymbol{J}}\begin{bmatrix}\dot{p}\\\ \dot{q}\\\ \dot{r}\\\
\end{bmatrix}=\begin{bmatrix}{L(T_{2}-T_{4})}\\\ {L(T_{3}-T_{1})}\\\
{M_{1}-M_{2}+M_{3}-M_{4}}\end{bmatrix}-\begin{bmatrix}p\\\ q\\\ r\\\
\end{bmatrix}\times{\boldsymbol{J}}\begin{bmatrix}p\\\ q\\\ r\\\
\end{bmatrix}$ (9)

From Eq. 8 and Eq. 9, the total control actuation inputs, defined as the
vecotr $\boldsymbol{U}=[u_{1},u_{2},u_{3},u_{4}]^{\top}$,

$\boldsymbol{U}=\begin{bmatrix}u_{1}\\\ u_{2}\\\ u_{3}\\\ u_{4}\\\
\end{bmatrix}=\begin{bmatrix}1&1&1&1\\\ 0&L&0&-L\\\ -L&0&L&0\\\
\gamma&-\gamma&\gamma&-\gamma\end{bmatrix}\begin{bmatrix}T_{1}\\\ T_{2}\\\
T_{3}\\\ T_{4}\\\ \end{bmatrix}$ (10)

where ${L}$ is the distance from the rotor axis of rotation to the quadrotor
center of mass and $\gamma=\frac{K_{M}}{K_{F}}$, where ${K_{m}}$ and ${K_{t}}$
are respectively the aerodynamic motor moment and force constants [42, Section
3.2].

###  2.4 Simulation

For the scope of this manuscript, we will partially focus on the simulation of
a 3 DOF quadrotor system, as shown in Figure 2 Quad. This design, also known
as a planar quadrotor, restricts the quadrotor’s motion to the $y-z$ plane
exclusively. At first, we wanted to use a simplified model to validate and
test SINDy algorithm to the problem in hand. Then, gradually we transition to
the entire 6-DOF analysis by studying this simplified subset derived from the
comprehensive 6-DOF scenario mentioned in Section 2.2. The simulation of the
system’s equations of motion will be done discretely using the well-known
Runge-Kutta fourth-order (RK4) numerical technique.

$\boldsymbol{x}_{k+1}=\boldsymbol{f}_{\text{RK4}}(\boldsymbol{x}_{k},\boldsymbol{u}_{k},\Delta
t)$ (11)

where $k$ is discrete time index, $\Delta t$ is the time step,
$\boldsymbol{u}$ is the control actuation. For the full picture formulation of
RK4 readers may refer to [43].

The 3 DOF equations of motion are

$\displaystyle\ddot{y}$ $\displaystyle=-\frac{u_{1}}{m}\sin{\phi}$ (12)
$\displaystyle\ddot{z}$ $\displaystyle=\frac{u_{1}}{m}\sin{\phi}-g$
$\displaystyle\ddot{\phi}$ $\displaystyle=\frac{u_{2}}{J_{x}}$

Figure 2: 3 DOF quadrotor with the inertial and the body reference frames.

or in a matrix format as

$\left[\begin{array}[]{c}\ddot{y}\\\ \ddot{z}\\\
\ddot{\phi}\end{array}\right]=\left[\begin{array}[]{c}0\\\ -g\\\
0\end{array}\right]+\left[\begin{array}[]{cc}-\frac{1}{m}\sin(\phi)&0\\\
\frac{1}{m}\cos(\phi)&0\\\
0&\frac{1}{I_{xx}}\end{array}\right]\left[\begin{array}[]{l}u_{1}\\\
u_{2}\end{array}\right]$ (13)

We consider the state vector
$\boldsymbol{x}=\begin{bmatrix}y,z,\phi,\dot{y},\dot{z},\dot{\phi}\end{bmatrix}^{\top}$.
Then, Equations 13 are linearized about the hovering position using the first
order Taylor series expansion approximation of the non-linear terms and
developed PD controller for the robot.

The linearized version of Eqs. 13

$\displaystyle\ddot{y}$ $\displaystyle=-g\phi$ (14) $\displaystyle\ddot{z}$
$\displaystyle=-g+\frac{u_{1}}{m}$ $\displaystyle\ddot{\phi}$
$\displaystyle=\frac{u_{2}}{I_{xx}}$

Introducing a generic state variable $\boldsymbol{\varrho}$ that can represent
either $\boldsymbol{y}$, $\boldsymbol{z}$, or $\boldsymbol{\phi}$ at a time.
For this generic state variable $\boldsymbol{\varrho}$ to be driven to a new
desired state, the commanded acceleration $\boldsymbol{\ddot{\varrho}}_{c}$ is
needed. For this we define the proportional and derivative error as

$\displaystyle e_{p}=\varrho_{d}-\varrho$ (15) $\displaystyle
e_{v}=\dot{\varrho}_{d}-\dot{\varrho}$

where the $\varrho_{d}$ is the desired value. Then we want both $e_{p}$ and
$e_{v}$ to satisfy,

$\left(\ddot{\varrho}_{d}-\ddot{\varrho}_{c}\right)+k_{p}e_{p}+k_{v}e_{v}=0$
(16)

in order to guarantee the error’s convergence under some values of $k_{p}$ and
$k_{v}$. So

$\displaystyle
u_{1}=mg+m\ddot{z}_{c}=m\left(\ddot{z}_{d}+k_{v,z}\left(\dot{z}_{d}-\dot{z}\right)+k_{p,z}\left(z_{d}-z\right)+g\right)$
(17) $\displaystyle
u_{2}=J_{x}\ddot{\phi}_{d}=J_{x}\left(\ddot{\phi}_{c}+k_{v,\phi}\left(\dot{\phi}_{c}-\dot{\phi}\right)+k_{p,\phi}\left(\phi_{c}-\phi\right)\right)$
$\displaystyle\phi_{c}=-\frac{\ddot{y}_{c}}{g}=-\frac{1}{g}\left(\ddot{y}_{d}+k_{v,y}\left(\dot{y}_{d}-\dot{y}\right)+k_{p,y}\left(y_{d}-y\right)\right)$

From this, we can simulate the 3 DOF quadrotor based on different desired
trajectory cases found in Table 1.

We have set the mass of the quadrotor to be $0.18$ Kg, the arm length is 0.086
m, and the quadrotor’s moment of inertia ($J_{x}$) is 0.00025 Kg.m2.

Case | Trajectory type | Initial states | Desired states
---|---|---|---
A | step | $x_{0}=[0,\;0,\;0,\;0,\;0,\;0]^{\top}$ | $x_{d}=[0.5,\;0.2,\;0,\;0,\;0,\;0]^{\top}$
B | Sine | $x_{0}=[0,\;0,\;0,\;0,\;0,\;0]^{\top}$ | $x_{d}=[4,\;0,\;0,\;0,\;0,\;0]^{\top}$
C | Diamond | $x_{0}=[0,\;1.8,\;0,\;0,\;0,\;0]^{\top}$ | $x_{d}=[0,\;0,\;0,\;0,\;0,\;0]^{\top}$
Table 1: Simulation cases for different trajectories.

The response of the three cases is summarized in Figures 3, 4, and 5. In the
simulation settings, one can see the acceptable tracking behaviour of the
quadrotor. The resulted data from will be used as a training and test data for
SINDy algorithm; more on this in section 3.1.

Figure 3: The reference step trajectory represent case A and the system
response Figure 4: The reference sinusoidal trajectory represent case B and
the system response Figure 5: The reference diamond-shaped trajectory
represent case C and the system response

##  3 Methodology

###  3.1 SINDy With Control Framework

Nonlinear dynamical systems can be represented as

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$ (18)

and we consider more general case when control takes place in the nonlinear
dynamics described as

$\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x},\boldsymbol{u})$ (19)

where $\boldsymbol{x}\in\mathbb{R}^{n}$ is the system states vector,
$\boldsymbol{u}\in\mathbb{R}^{q}$ is the control vector, and $\boldsymbol{f}$
is a nonlinear such that
$\boldsymbol{f}:\mathbb{R}^{n}\times\mathbb{R}^{q}\longrightarrow\mathbb{R}^{n}$.

SINDy with control technique utilizes sparse regression to identify a minimal
set of active terms from a candidate library
$\boldsymbol{\Pi(\boldsymbol{x},\boldsymbol{u})}$ of linear and nonlinear
terms in the system states $\boldsymbol{x}$ and actuation variables
$\boldsymbol{u}$ that approximate the underlying nonlinear dynamics
represented by $\boldsymbol{f}$. On the premise that a considerable number of
systems manifest relatively sparse active terms in their dynamics. To
construct $\boldsymbol{\Pi}$ we measure $m$ snapshots of the state
$\boldsymbol{x}$ and actuation $\boldsymbol{u}$ by experiment or from
numerical simulation.

$\boldsymbol{X}=\begin{bmatrix}\boldsymbol{x}^{\top}(t_{1})\\\
\boldsymbol{x}^{\top}(t_{2})\\\ \boldsymbol{x}^{\top}(t_{3})\\\ \vdots\\\
\boldsymbol{x}^{\top}(t_{m})\end{bmatrix}=\begin{bmatrix}x_{1}(t_{1})&x_{2}(t_{1})&\dots&x_{n}(t_{1})\\\
x_{1}(t_{2})&x_{2}(t_{2})&\dots&x_{n}(t_{2})\\\
x_{1}(t_{3})&x_{2}(t_{3})&\dots&x_{n}(t_{3})\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{1}(t_{m})&x_{2}(t_{m})&\dots&x_{n}(t_{m})\end{bmatrix},\quad\boldsymbol{U}=\begin{bmatrix}\boldsymbol{u}^{\top}(t_{1})\\\
\boldsymbol{u}^{\top}(t_{2})\\\ \boldsymbol{u}^{\top}(t_{3})\\\ \vdots\\\
\boldsymbol{u}^{\top}(t_{m})\end{bmatrix}=\begin{bmatrix}u_{1}(t_{1})&u_{2}(t_{1})&\dots&u_{q}(t_{1})\\\
u_{1}(t_{2})&u_{2}(t_{2})&\dots&u_{q}(t_{2})\\\
u_{1}(t_{3})&u_{2}(t_{3})&\dots&u_{q}(t_{3})\\\ \vdots&\vdots&\ddots&\vdots\\\
u_{1}(t_{m})&u_{2}(t_{m})&\dots&u_{q}(t_{m})\end{bmatrix}$ (20)

The matrix $\boldsymbol{\Pi}$ can be reconstructed now by

$\boldsymbol{\Pi}(\boldsymbol{X,U})=\begin{bmatrix}|&|&|&|&|&|&&|&|&\\\
\boldsymbol{1}&\boldsymbol{X}&\boldsymbol{U}&\boldsymbol{X\otimes
X}&\boldsymbol{X\otimes U}&\boldsymbol{U\otimes
U}&\dots&\cos(\boldsymbol{X})&\cos(\boldsymbol{X\otimes X})&\dots\\\
|&|&|&|&|&|&&|&|&\end{bmatrix}$ (21)

The effectiveness of the candidate term library is important to the SINDy with
control algorithm. A basic strategy starts with a simple option, like
polynomials, and gradually increases the complexity of the library by
incorporating additional terms [3]. After evaluating the library, the states
derivatives also should be evaluated. Here we evaluated the derivatives by
using numerical approximation. Specifically, we used the finite difference
method.

The system in Eq. (19) can be rewritten as

$\dot{\boldsymbol{X}}=\boldsymbol{\Pi}(\boldsymbol{X},\boldsymbol{U})\boldsymbol{\Omega}$
(22)

where $\boldsymbol{\Lambda}$ is a coefficient matrix that is almost sparse.
This matrix tries to activate the fewest active terms in the candidate matrix
$\boldsymbol{\Pi}$ that results in the best model fit:

$\omega_{j}=\operatorname*{argmin}_{\tilde{\omega}_{j}}||\dot{\boldsymbol{X}}-\boldsymbol{\Pi}(\boldsymbol{X,U})\tilde{\omega}_{j}||_{2}^{2}+\lambda||\tilde{\omega}_{j}||_{1}$
(23)

This problem can be solved using variety of regression techniques. Including
but not limited to, LASSO [27], sequential thresholded least-squares [7], and
Sparse Relaxed Regularized Regression (SR3) [44]. We tried different
optimizers and have found SR3 is the best for our case.

####  3.1.1 Hyperparameter Tuning

The hyperparameter $\lambda$ is critical in identifying the most sparse
dynamics. Typically, $\lambda$ takes values between 0 to 1. So, we scanned
over the lambda domain with step size of 0.05 and evaluated the model using
test data to cross validate utilizing the root mean squared error (RMSE). We
found the best value of $\lambda$ to be equal to 0.45.

We assume the generic prediction of a generic state vector
$\boldsymbol{\varrho}$ is denoted as $\hat{\boldsymbol{\varrho}}$ and the RMSE
of $\boldsymbol{\varrho}$ is denoted as $\tilde{\varrho}$

$\tilde{\varrho}=\sqrt{\frac{\sum_{i=1}^{N}(\varrho_{i}-\hat{\varrho}_{i})^{2}}{N}}$
(24)

The chosen lambda corresponds to the RMSE of the states of the planar quad
case.

State | ${y}$ | ${z}$ | ${\phi}$ | $\dot{{y}}$ | $\dot{{z}}$ | $\dot{{\phi}}$
---|---|---|---|---|---|---
RMSE $\times 10^{-3}$ | 0.0088 | 0.0152 | 0.0227 | 0.0294 | 0.0070 | 0.1379
Table 2: The RMSE for the chosen $\lambda=0.45$ for diamond shaped trajectory.

####  3.1.2 Candidate functions

In the present study, we exploit our comprehensive understanding of the
physical system in hand. Specifically, we propose that the presence of
nonlinearity in the system can be expressed through polynomials and Fourier
basis functions, such as those involving sine and cosine. Through mathematical
analysis of the system, we have discerned that from the entire state space of
the system, only the Euler angles can be represented as Fourier basis
functions, and the rest can be characterized by polynomials.

###  3.2 Training

The numerical data from the 3 DOF planar quadrotor simulation was utilized. We
first choose Case C because we thought it would allow the quadrotor to span
the entire state space. Consequently, it will be the ideal case to start with
for training. It was subsequently demonstrated by the other cases that the
algorithm is generic regardless of trajectory. As long as the trajectory gave
sufficient data and spanned the state space appropriately, the algorithm
successfully captured the complete dynamics. If, on the other hand, the
quadrotor followed a signal that led it to move only in a straight line along
the $y$ or $z$ axis, the algorithm may fail to distinguish the unseen states,
resulting in an incomplete identification of the dynamics. We used snapshots
with $\Delta t=0.05$ with 1000 snapshots in time. We differentiated the states
using finite difference numerical scheme of order one.

##  4 Results and Discussion

###  4.1 Discovered Model

The model is trained with the data extracted from case C as discussed earlier
and it came out with the following dynamics

$\displaystyle\dot{y}$ $\displaystyle=\dot{y}$ (25) $\displaystyle\dot{z}$
$\displaystyle=\dot{z}$ $\displaystyle\dot{\phi}$
$\displaystyle=0.993\dot{\phi}$ $\displaystyle\ddot{y}$
$\displaystyle=-5.549u_{1}\sin(\phi)$ $\displaystyle\ddot{z}$
$\displaystyle=-9.811+5.556u_{1}\cos(\phi)$ $\displaystyle\ddot{\phi}$
$\displaystyle=4000.000u_{2}$

That nearly matches the original derived mathematical model from the first
principles. Table 3 shows a comparison between the discovered dynamics by
SINDy and the gound truth mathematical model.

Table 3: Comparison between the discovered dynamics and ground truth mathematical model. States | SINDy | Mathematical Model
---|---|---
$\dot{y}$ | $1.000\dot{y}$ | $1.000\;\dot{y}$
$\dot{z}$ | $1.000\;\dot{z}$ | $1.000\;\dot{z}$
$\dot{\phi}$ | $0.993\;\dot{\phi}$ | 1.000 $\dot{\phi}$
$\ddot{y}$ | $-5.549\;u_{1}\sin(\phi)$ | $-5.5556\;u_{1}\;\sin(\phi)$
$\ddot{z}$ | $-9.811+5.556\;u_{1}\cos(\phi)$ | $-9.81+5.556\;u_{1}\cos(\phi)$
$\ddot{\phi}$ | $4000.000\;u_{2}$ | $4000.000\;u_{2}$

However, we also used the other cases to train the model and it resulted in
models that compare to the original one as we discussed in section 3.2.

###  4.2 Testing

####  4.2.1 Case A

Here we simulate the discovered dynamics using the step trajectory as desired
trajectory. We compare the system behaviour over time for both the discovered
SINDy dynamics and the ground truth. Figure 6(b) shows the absolute error
between the predicted states $\hat{\boldsymbol{\varrho}}$ and the original
states $\boldsymbol{\varrho}$. The results show that the absolute error
approaches zero, giving strong validation for the identified model. This shows
that the SINDy algorithm accurately captures the underlying dynamics of the
system, closely matching the ground truth dynamics. The relatively small error
between anticipated and original states confirms the discovered model’s
efficacy and reliability in capturing the main characteristics of quadrotor
dynamics. The close agreement between the identified dynamics and the ground
truth confirms the SINDy algorithm’s utility as a good tool for extracting
mathematical models from data.

(a) The reference step trajectory represent case A and the system response

(b) The absolute error between the predicted states
$\hat{\boldsymbol{\varrho}}$ and the original states $\boldsymbol{\varrho}$

Figure 6: SINDy comparison with the mathematical model.

####  4.2.2 Case B

As in the previous section, we simulate the discovered dynamics using the
sinusoidal trajectory as desired trajectory. We compare the system behaviour
over time for both the discovered SINDy dynamics and the ground truth. Figure
7(b) shows the absolute error between the predicted states
$\hat{\boldsymbol{\varrho}}$ and the original states $\boldsymbol{\varrho}$.
We can say that the error nearly equal zero. This validates the discovered
model.

(a) The reference sinusoidal trajectory represent case A and the system
response

(b) The absolute error between the predicted states
$\hat{\boldsymbol{\varrho}}$ and the original states $\boldsymbol{\varrho}$

Figure 7: SINDy comparison with the mathematical model.

##  5 Conclusion

In this study, we have employed the SINDy algorithm to extract a mathematical
model from simulation data of a quadrotor. Initially, the quadrotor dynamics
were modeled in the 6 DOF configuration. Subsequently, for the purpose of
simplification and model validation, the system was constrained to a 3 DOF
representation within a 2-dimensional plane. To assess the effectiveness of
the SINDy algorithm, three distinct simulation cases, as depicted in Figure
LABEL:fig:response, were considered. A systematic exploration of the
hyperparameter $\lambda$ space was conducted. Through this thorough analysis,
we successfully identified the optimal model with an associated $\lambda$
value of 0.45, as demonstrated in Table 2, based on the RMSE metric. Table 3
and Figures 6(b), 7(b) show that the algorithm captured almost the same
dynamics as the quadrotor ground truth mathematical model.

We also demonstrated the SINDy powerfulness in obtaining solid mathematical
models from data, which will aid in the advancement of modeling in the field
of aerospace applications, particularly UAVs. Furthermore, in a future version
of this research, we will attempt to analyze the algorithm’s resilience in the
face of noisy measurements in order to assure its robustness in real-world
circumstances. Understanding how measurement errors affect the identified
model’s accuracy and reliability. Given the occurrence of mathematically
unmodeled disturbances and uncertainties in real-world scenarios, future
research should concentrate on building discrepancy models capable of properly
capturing and accounting for these aspects. Incorporating these models with
the SINDy algorithm can improve its ability to deal with disturbances,
resulting in more robust and accurate predictions.

Furthermore, we intend to extend the SINDy algorithm’s applicability to the
full 6 DOF quadrotor model, encompassing all elements of its dynamic behavior.
A more thorough understanding of the quadrotor’s dynamics can be obtained by
including the entire model, including translational and rotational motion.
There will also be an in-depth examination of the optimizers used in the SINDy
algorithm. Examining various optimization strategies and their impact on the
quality of the identified model can result in improvements in convergence
speed, accuracy, and the capacity to handle complicated datasets successfully.

## References

  * [1] D. K. Arrowsmith, C. M. Place, C. H. Place, and Arrowsmith/Place.  An Introduction to Dynamical Systems.  Cambridge University Press, July 1990.  Google-Books-ID: HOU6q6wHCncC.
  * [2] Steven L. Brunton and J. Nathan Kutz.  Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control.  Cambridge University Press, Cambridge, 2019.
  * [3] Eurika Kaiser, J Nathan Kutz, and Steven L Brunton.  Sparse identification of nonlinear dynamics for model predictive control in the low-data limit.  Proceedings of the Royal Society A, 474(2219):20180335, 2018.
  * [4] Clarence W Rowley, Tim Colonius, and Richard M Murray.  Model reduction for compressible flows using pod and galerkin projection.  Physica D: Nonlinear Phenomena, 189(1-2):115–129, 2004.
  * [5] Clarence W Rowley and Scott TM Dawson.  Model reduction for flow analysis and control.  Annual Review of Fluid Mechanics, 49:387–417, 2017.
  * [6] Mikhail Hayhoe, Francisco Barreras, and Victor M Preciado.  Data-driven control of the covid-19 outbreak via non-pharmaceutical interventions: A geometric programming approach.  arXiv preprint arXiv:2011.01392, 2020.
  * [7] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, April 2016.  Publisher: Proceedings of the National Academy of Sciences.
  * [8] Josh Bongard and Hod Lipson.  Automated reverse engineering of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 104(24):9943–9948, 2007.
  * [9] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  science, 324(5923):81–85, 2009.
  * [10] John R Koza.  Genetic programming as a means for programming computers by natural selection.  Statistics and computing, 4:87–112, 1994.
  * [11] Peter Schmid and Joern Sesterhenn.  Dynamic mode decomposition of numerical and experimental data.  In APS Division of Fluid Dynamics Meeting Abstracts, volume 61, pages MR–007, August 2008.
  * [12] B. O. Koopman.  Hamiltonian Systems and Transformation in Hilbert Space.  Proceedings of the National Academy of Sciences, 17(5):315–318, May 1931.  Publisher: Proceedings of the National Academy of Sciences.
  * [13] Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley.  A Data–Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition.  J Nonlinear Sci, 25(6):1307–1346, December 2015.
  * [14] Matthew O. Williams, Clarence W. Rowley, and Ioannis G. Kevrekidis.  A Kernel-Based Approach to Data-Driven Koopman Spectral Analysis, July 2015.  arXiv:1411.2260 [math].
  * [15] Peter J Baddoo, Benjamin Herrmann, Beverley J McKeon, and Steven L Brunton.  Kernel learning for robust dynamic mode decomposition: linear and nonlinear disambiguation optimization.  Proceedings of the Royal Society A, 478(2260):20210830, 2022.
  * [16] Steven L. Brunton, Bingni W. Brunton, Joshua L. Proctor, Eurika Kaiser, and J. Nathan Kutz.  Chaos as an intermittently forced linear system.  Nat Commun, 8(1):19, May 2017.  Number: 1 Publisher: Nature Publishing Group.
  * [17] Hassan Arbabi and Igor Mezić.  Ergodic theory, Dynamic Mode Decomposition and Computation of Spectral Properties of the Koopman operator.  SIAM J. Appl. Dyn. Syst., 16(4):2096–2126, January 2017.  arXiv:1611.06664 [math].
  * [18] Joshua L. Proctor, Steven L. Brunton, and J. Nathan Kutz.  Dynamic Mode Decomposition with Control.  SIAM J. Appl. Dyn. Syst., 15(1):142–161, January 2016.  Publisher: Society for Industrial and Applied Mathematics.
  * [19] Jonathan H. Tu, Clarence W. Rowley, Dirk M. Luchtenburg, Steven L. Brunton, and J. Nathan Kutz.  On Dynamic Mode Decomposition: Theory and Applications.  Journal of Computational Dynamics, 1(2):391–421, 2014.  arXiv:1312.0041 [physics].
  * [20] Jonathan H Tu.  Dynamic mode decomposition: Theory and applications.  PhD thesis, Princeton University, 2013.
  * [21] Kathleen Champion.  From data to dynamics: discovering governing equations from data.  PhD thesis, 2019.
  * [22] Enoch Yeung, Soumya Kundu, and Nathan Hodas.  Learning Deep Neural Network Representations for Koopman Operators of Nonlinear Dynamical Systems, November 2017.  arXiv:1708.06850 [cs, math].
  * [23] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi.  Learning Koopman Invariant Subspaces for Dynamic Mode Decomposition, January 2018.  arXiv:1710.04340 [cs, math, stat].
  * [24] E.J. Candes, J. Romberg, and T. Tao.  Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information.  IEEE Transactions on Information Theory, 52(2):489–509, February 2006.  Conference Name: IEEE Transactions on Information Theory.
  * [25] Emmanuel J. Candès, Justin K. Romberg, and Terence Tao.  Stable signal recovery from incomplete and inaccurate measurements.  Communications on Pure and Applied Mathematics, 59(8):1207–1223, 2006.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.20124.
  * [26] David Donoho.  Compressed sensing.  IEEE Transactions on Information Theory, 52(4), April 2006.
  * [27] Robert Tibshirani.  Regression Shrinkage and Selection Via the Lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1996.tb02080.x.
  * [28] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani.  An Introduction to Statistical Learning: with Applications in R.  Springer Texts in Statistics. Springer US, New York, NY, 2021.
  * [29] Patrick Gelß, Stefan Klus, Jens Eisert, and Christof Schütte.  Multidimensional approximation of nonlinear dynamical systems.  Journal of Computational and Nonlinear Dynamics, 14(6), 2019.
  * [30] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Data-driven discovery of partial differential equations.  Science advances, 3(4):e1602614, 2017.
  * [31] Hayden Schaeffer.  Learning partial differential equations via data discovery and sparse optimization.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2197):20160446, 2017.
  * [32] Hayden Schaeffer and Scott G. McCalla.  Sparse model selection via integral terms.  Physical Review E, 96(2):undefined–undefined, 2017.  Number: 2.
  * [33] Patrick AK Reinbold, Daniel R Gurevich, and Roman O Grigoriev.  Using noisy or incomplete data to discover models of spatiotemporal dynamics.  Physical Review E, 101(1):010203, 2020.
  * [34] Yu Xin Jiang, Xiong Xiong, Shuo Zhang, Jia Xiang Wang, Jia Chun Li, and Lin Du.  Modeling and prediction of the transmission dynamics of COVID-19 based on the SINDy-LM method.  Nonlinear Dynamics, 105(3):2775–2794, 2021.  Number: 3.
  * [35] Gregor Thiele, Arne Fey, David Sommer, and Jorg Kruger.  System identification of a hysteresis-controlled pump system using SINDy.  2020 24th International Conference on System Theory, Control and Computing, ICSTCC 2020 - Proceedings, pages 457–464, 2020.
  * [36] Bhavana Bhadriraju, Mohammed Saad Faizan Bangi, Abhinav Narasingam, and Joseph Sang Il Kwon.  Operable adaptive sparse identification of systems: Application to chemical processes.  AIChE Journal, 66(11):undefined–undefined, 2020.  Number: 11.
  * [37] Dipankar Bhattacharya, Leo K. Cheng, and Weiliang Xu.  Sparse Machine Learning Discovery of Dynamic Differential Equation of an Esophageal Swallowing Robot.  IEEE Transactions on Industrial Electronics, 67(6):4711–4720, 2020.  Number: 6.
  * [38] Mariia Sorokina, Stylianos Sygletos, and Sergei Turitsyn.  Sparse Identification for Nonlinear Optical Communication Systems: SINO Method.  Opt. Express, 24(26):30433, December 2016.  arXiv:1701.01650 [physics].
  * [39] Tayyab Manzoor, Hailong Pei, Zhongqi Sun, and Zihuan Cheng.  Model predictive control technique for ducted fan aerial vehicles using physics-informed machine learning.  Drones, 7(1):4, 2022.
  * [40] Markus Quade, Markus Abel, J Nathan Kutz, and Steven L Brunton.  Sparse identification of nonlinear dynamics for rapid model recovery.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(6):063116, 2018.
  * [41] Kadierdan Kaheman, Eurika Kaiser, Benjamin Strom, J Nathan Kutz, and Steven L Brunton.  Learning discrepancy models from experimental data.  arXiv preprint arXiv:1909.08574, 2019.
  * [42] HTMN ElKholy.  Dynamic modeling and control of a quadrotor using linear and nonlinear approaches, 2014.  Master thesis.
  * [43] Joan Sola.  Quaternion kinematics for the error-state kalman filter.  arXiv preprint arXiv:1711.02508, 2017.
  * [44] Peng Zheng, Travis Askham, Steven L Brunton, J Nathan Kutz, and Aleksandr Y Aravkin.  A unified framework for sparse relaxed regularized regression: Sr3.  IEEE Access, 7:1404–1423, 2018.


# The third law of thermodynamics and black holes

H<EMAIL_ADDRESS>A. H.
<EMAIL_ADDRESS>Iarley P.
<EMAIL_ADDRESS>J. P. Morais<EMAIL_ADDRESS>U. K<EMAIL_ADDRESS>A. Sayahian Jahromi5 1 Research
Institute for Astronomy and Astrophysics of Maragha (RIAAM), University of
Maragheh, P.O. Box 55136-553, Maragheh, Iran
2 Department of Chemistry and Physics, Federal University of Paraíba, Rodovia
BR 079 - Km 12, 58397-000 Areia-PB, Brazil
3 Instituto de Física, Universidade Federal do Rio de Janeiro, 21.941-972 -
Rio de Janeiro-RJ-Brazil
4 Department of Mathematics, Institute of Applied Sciences and Humanities, GLA
University, Mathura-281406, Uttar Pradesh, India
5 Zarghan Branch, Islamic Azad University, Zarghan, Iran

###### Abstract

Working in the framework of generalized statistics, the problem of
availability of the third law of thermodynamics in the black hole physics is
studied by focusing on Schwarzschild black hole which easily and clearly
exposes the violation of this law in the common approach based on Bekenstein
entropy. Additionally, it is addressed that some inconsistencies between the
predictions of quantum field theory and those of thermodynamics on the black
hole temperature may be reconciled by using the thermodynamics laws in order
to broaden energy definition. It claims that thermodynamics should be employed
as a powerful tool in looking for more comprehensive energy definitions in
high energy physics, still mysterious.

….

###### pacs:

….

##  I Introduction

The third law of thermodynamics states that the entropy of a system should
approach a constant value ($C$) at absolute zero temperature, or equally
$S(T\rightarrow 0)\rightarrow C$. Bekenstein entropy ($S_{B}$) is proportional
to the horizon area $(A=4\pi r_{h}^{2})$, and correspondingly, for a
Schwarzschild black hole of mass $M(\equiv E)$ and Hawking temperature $T_{H}$
for which $r_{h}=2M$ and $T=\frac{\partial M}{\partial S}=\frac{1}{8\pi
M}=T_{H}$, we have

$\displaystyle S_{B}=\frac{A}{4}=4\pi M^{2}=\frac{1}{16\pi T_{H}^{2}},$ (1)

in the units of $c=\hbar=k_{B}=G=1$, where $k_{B}$ denotes the Boltzmann
constant. It clearly indicates that the third law of thermodynamics is not
satisfied or briefly, $T_{H}\rightarrow 0(\parallel
M\rightarrow\infty)\Rightarrow S_{B}\rightarrow\infty$. Although we considered
Schwarzschild black hole, the behavior of
$S_{B}\big{(}T(M\rightarrow\infty)\rightarrow 0\big{)}\rightarrow\infty$ is
common in other black holes such as Kerr-Newman and Reissner-Nordström metrics
1 ; 01 ; 2 ; 3 ; 4 .

$S_{B}$ is non-extensive, a trait which reminds the generalized statistics
such as Tsallis statistics revT ; masi  including entropies which are not
extensive. Indeed, as gravitational systems include long-range interaction
(gravity), it is proposed that the Boltzmann entropy (leading to Eq. (1))
should be replaced with generalized entropies which leads to substantial
consequences in both gravitational and cosmological setups (see for example
Refs. tsallis ; refgerg ; gerg ; non13 ; nonK ; EPJC ; KHDE ; sadeghi ; mesri2
their references and citations). It also seems that there is a connection
between deviation from the Boltzmann statistics and the quantum features of
gravity epl ; homa ; mesri ; barrow ; mesri2 , and consequently, the
generalized and loop quantum gravity entropies can be classified as the
subclasses of a general entropy gen .

Our first aim is to study the possibility of satisfying the third law by
employing some new entropies proposed that provide respectable solutions in
cosmological and gravitational setups. Hereof and in subsequent sections, we
respectively study the problem by considering Tsallis and Cirto entropy,
Tsallis entropy and Kaniadakis entropy. Throughout the survey, we also address
a thermodynamic energy definition for a black hole of mass $M$, corresponding
to each entropy, which admits Hawking temperature. The black hole remnant and
its decay time in mentioned entropy formalisms have been studied in fifth
section. The last section is also devoted to summary.

##  II Tsallis and Cirto entropy and the third law of thermodynamics

Motivated by the non-extensivity of $S_{B}$, and also the long-range nature of
gravity, Tsallis and Cirto tsallis  introduce a new entropy for black holes as

$\displaystyle S_{T}=\gamma A^{\delta},$ (2)

where $\gamma$ and $\delta$ are two unknown free constants evaluated by other
parts of physics or observations. It is also useful to note here that this
form of entropy is also supported in the framework of quantum gravity barrow .
It means that two different approaches and motivations lead to the same result
which increases the credit of this new proposal for the black hole entropy.
Moreover, the equality of results can be considered as the sign of a deep
connection between non-extensivity and quantum gravity helping us build a
relation between their parameters, a result noted in Refs. epl ; homa .
Considering $r_{h}=2M$ and $T=\frac{\partial M}{\partial
S_{T}}=\frac{1}{2\delta\gamma(16\pi)^{\delta}M^{2\delta-1}}$, one easily finds
that the third law of thermodynamics is met whenever $0<\delta<\frac{1}{2}$,
and in summary, $S_{T}\rightarrow 0\parallel M\rightarrow 0\parallel
T\rightarrow 0$.

Now, let us employ Hawking temperature ($T_{H}=\frac{1}{8\pi M}$) instead of
$T=\frac{\partial M}{\partial S_{T}}$ which leads to $S_{T}\propto
T^{-2\delta}$ meaning that the third law is fulfilled only if $\delta<0$ and
for this case, we briefly have $M\rightarrow\infty\parallel S_{T},T\rightarrow
0$. In Tsallis statistics, an intrinsic temperature discrepancy between real
temperature and the temperature obtained by thermodynamic relation may emerge
depending on the expectation values definition (averaging methods) used in
obtaining quantities such as energy kol , and only, having in hand the system
temperature, one can decide on true temperature kol . Therefore, the above
obtained temperature discrepancy may be more understandable by bearing in mind
the intrinsic temperature discrepancy of Tsallis statistics. Consequently,
since $\delta$ is a free parameter estimated from observations revT ; masi  or
probably other parts of physics epl ; homa , we cannot go further in choosing
one of these temperatures, and thus the corresponding thermodynamics, unless
we enclasp detailed observations, data and info on black holes.

Of course, a way to reconcile the above inconsistency between temperatures is
to redefine energy. In both cases above, we assumed $E=M$, while if we assume
$T=T_{H}=\frac{\partial E}{\partial S}$, and use Eq. (2), then we reach

$\displaystyle E_{T}=\int_{0}^{M}\frac{1}{8\pi m}\frac{\partial
S_{T}}{\partial m}dm,$ (3)

finally leading to

$\displaystyle
E_{T}=4\gamma\delta\frac{(4\pi)^{\delta-1}}{2\delta-1}M^{2\delta-1},$ (4)

as the energy of a Schwarzschild black hole of mass $M$ in Tsallis formalism
which recovers $E=M$ by inserting $\delta=1$ and $\gamma=\frac{1}{4}$ (the
Bekenstein limit). It is also obvious that $E_{T}$ is positive if
$\delta>\frac{1}{2}$ or $\delta<0$. For the $\delta>\frac{1}{2}$ case, the
third law is not satisfied and we face a situation like to what we obtained in
the case of Bekenstein entropy (i. e, $T_{H}\rightarrow 0\parallel
S_{B},M,E_{T}\rightarrow\infty$). The third law is met for $\delta<0$ and in
parallel $E_{T}\rightarrow 0$ or briefly, $E_{T},S_{T},T\rightarrow 0\parallel
M\rightarrow\infty$.

Finally, we think that since the Hawking temperature is also supported by
other parts of physics such as quantum field theory in curved spacetime and so
on haw ; bek ; thes , and as there is not any common agreement on the energy
definition in high energy physics energy ; energy1 , it is probably reasonable
to rely on this approach compared to the two previously mentioned cases. The
latter means that thermodynamics may be used to find a more proper energy
definition in high energy physics.

In summary, we obtained $3$ cases recaped as

$\displaystyle\\!\\!\\!\\!i)\ E=M,\ T=\frac{\partial E}{\partial
S_{T}}=\frac{1}{2\delta\gamma(16\pi)^{\delta}M^{2\delta-1}}\neq T_{H},$
$\displaystyle\\!\\!\\!\\!ii)\ E=M,\ T=T_{H}=\frac{1}{8\pi
M}\neq\frac{\partial E}{\partial S_{T}},$ (5)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!iii)\
E_{T}=4\gamma\delta\frac{(4\pi)^{\delta-1}}{2\delta-1}M^{2\delta-1},\
T=T_{H}=\frac{1}{8\pi M}=\frac{\partial E_{T}}{\partial S_{T}},$

where the third law is satisfied for the first case when
$0<\delta<\frac{1}{2}$, and for the remaining cases when $\delta<0$.

##  III Tsallis entropy and the third law

Recently, focusing on the relation between Tsallis and Boltzmann statistics, a
new entropy has been derived for black holes as KHDE

$\displaystyle S_{q}=\frac{1}{1-q}[\exp\big{(}(1-q)S_{B}\big{)}-1]$
$\displaystyle=\frac{2\exp(\frac{(1-q)S_{B}}{2})}{(1-q)}\sinh\left(\frac{(1-q)S_{B}}{2}\right),$
(6)

where $q$ is a free unknown parameter and this result is also confirmed by
calculating the Tsallis entropy content of black holes in the framework of
quantum gravity mesri ; KHDE . Following the recipe of previous section, if we
assume $E=M$ then we have

$\displaystyle T=\frac{\partial E}{\partial S_{q}}=T_{H}\exp\big{(}(q-1)4\pi
M^{2}\big{)},$ (7)

which recovers $T_{H}$ whenever $q=1$ (the Bekenstein limit of (III) KHDE ).
Briefly, $T\rightarrow 0$ only if $q<1$, and in this manner
$M,S_{q}\rightarrow\infty$ meaning that the third law is not satisfied. On the
other hand, if we assume $T=T_{H}$ (case $ii$), then we see that the third law
is satisfied only if $q>1$, or briefly $M\rightarrow\infty\Rightarrow
S_{q}(T\rightarrow 0)\rightarrow 0$. For the case $iii$, where
$T=T_{H}=\frac{1}{8\pi M}=\frac{\partial E_{q}}{\partial S_{q}}$, we reach

$\displaystyle E_{q}=\int_{0}^{M}\exp\big{(}(1-q)4\pi m^{2}\big{)}dm,$ (8)

as the energy content of a black hole of mass $M$. Clearly, the third law is
again satisfied only if $q>1$ and moreover, $E=\int_{0}^{M}dm=M$ is recovered
at the limit of $q\rightarrow 1$. The above integration can be performed with
the solution given as

$\displaystyle E_{q}=\frac{{\rm
erf}\left(2\sqrt{\pi}M\sqrt{q-1}\right)}{4\sqrt{q-1}},$ (9)

in which $\rm erf(x)$ denotes the error function er .

##  IV Kaniadakis entropy and the third law

Kaniadakis entropy of a black hole is also reported as KHDE

$\displaystyle S_{\kappa}=\frac{1}{\kappa}\sinh\big{(}\kappa S_{B}\big{)},$
(10)

where $\kappa$ is an unknown parameter evaluated by observations and probably,
the other parts of physics KHDE . Here, simple calculations lead to

$\displaystyle T=\frac{\partial E}{\partial
S_{\kappa}}=\frac{T_{H}}{\cosh(\kappa S_{B})},$ (11)

for the $i$-th case. The result indicates that, independent of the value of
$\kappa$, the third law is not satisfied
($S_{\kappa}\rightarrow\infty\parallel M\rightarrow\infty\parallel
T\rightarrow 0$). For the second case ($T=T_{H},E=M$), we can write

$\displaystyle S_{\kappa}=\frac{1}{\kappa}\sinh\left(\frac{\kappa}{16\pi
T^{2}}\right),$ (12)

which shows the third law is satisfied only if $\kappa<0$. If $\kappa<0$, then
the third law is also met for case $iii$, where $T=T_{H}$ and energy content
of black hole is obtainable by using

$\displaystyle E_{\kappa}=\int_{0}^{M}\cosh(\kappa 4\pi m^{2})dm,$ (13)

which recovers the $E=\int_{0}^{M}dm=M$ results at the limit of $\kappa=0$.
The solution to the above integral is also given as

$\displaystyle E_{\kappa}=\frac{1}{8\sqrt{\kappa}}\left[{\rm
erf}\left(2\sqrt{\kappa\pi}M\right)+{\rm
erfi}\left(2\sqrt{\kappa\pi}M\right)\right],$ (14)

where

${\rm erfi}(x)=-i{\rm erf}(ix).$ (15)

In Fig. (1), $E_{\kappa}$ and $E_{q}$ are plotted for different values of $q$
and $\kappa$, the $E=M$ case has also been depicted to have a comparison. It
is worthwhile to mention that, as it is obvious from this figure, there is an
asymptote for $E_{q}$ when $q>1$ as $E_{q}(M\gg
1)\rightarrow\frac{1}{4\sqrt{q-1}}$.

Figure 1: Upper panel: The behavior of energy content of the Schwarzschild
black hole as a function of its mass for different values of $q$ parameter.
Lower panel: The behavior of energy content of the Schwarzschild black hole as
a function of its mass for different values of $\kappa$ parameter. The dashed
curve also shows the Schwarzschild black hole of $E=M$.

##  V Black body radiation and black hole evaporation

In the framework of common statistical mechanics based on Gibbs entropy, the
black hole evaporation is described by Stefan-Boltzmann (SB) formula cav ; ali
. As it is apparent from Eq. (1), we have $S_{B},M\rightarrow\infty$ when
$T\rightarrow 0$ meaning that we face a catastrophe ali . This result can be
summarized in the language of SB law as

$\displaystyle\frac{dM}{dt}=-A\sigma T^{4},$ (16)

where the ordinary energy definition ($E=M$) is used and
$\sigma(=\frac{\pi^{2}}{60}$ in our units kanbb1 ; kelly ) denotes the SB
constant ali . Here, minus sign is also due to the fact that Eq. (16) explains
the amount of energy that system loses it. It clearly shows that the decay
rate ($\frac{dM}{dt}$) diverges as $M$ ($T$) approaches zero (infinity) ali .
Another consequence of this law is our ability to find the decay time $\tau$
as

$\displaystyle\tau=-\frac{1}{\sigma}\int_{M}^{0}\frac{dm}{AT^{4}}=\frac{(8\pi)^{3}}{2\sigma}\frac{M^{3}}{3},$
(17)

and therefore $\tau\sim M^{3}$ for a Schwarzschild black hole of temperature
$T_{H}$. Eq. (17) also indicates that $\tau\rightarrow\infty$ when
$M\rightarrow\infty$ while $T,\frac{dM}{dt}\rightarrow 0$ ali .

###  V.1 Tsallis and Cirto (TC) black hole

Black body radiation has not been studied in the formalism of Eq. (2), but as
we mentioned, this entropy is also confirmed by quantum gravity barrow . The
latter allows us to use the black body radiation formula in the framework of
quantum gravity to study the black holes meeting Eq. (2). The modifications of
quantum gravity on the black body spectrum have recently been studied by some
authors in different bases nozar ; hus ; cqg ; lobo . A common feature of
quantum gravity scenarios is generalized uncertainty principle (GUP) that
motivates us to focus on its modification on black body radiation. Indeed,
such a modification is also obtainable by using other quantum gravity
scenarios hus ; cqg .

In this regard, it is proposed that GUP may modify the black body spectrum as
cqg

$\displaystyle\frac{dE}{dt}=-A\sigma\big{[}T^{4}+\frac{15}{9}\alpha
T^{6}\big{]},$ (18)

where $\alpha$ is called the GUP parameter originated from the quantum aspects
of gravity, and there is an ongoing debate on its value agha ; cqg . It is
also obvious that Eq. (16) is recovered whenever $\alpha\rightarrow 0$ and
$E=M$. Calculations for three obtained cases lead to

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{i}^{TC}$ $\displaystyle=$
$\displaystyle-\frac{(2\delta\gamma(16\pi)^{\delta})^{4}}{16\pi\sigma}\int_{M}^{0}\frac{dm}{m^{6-8\delta}+\frac{15}{(6\delta\gamma(16\pi)^{\delta})^{2}}\alpha
m^{8-12\delta}},$ $\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{ii}^{TC}$
$\displaystyle=$
$\displaystyle-\frac{1}{18\sigma}\Big{[}1536\pi^{3}M^{3}-120\pi M\alpha$
$\displaystyle+$ $\displaystyle
5\sqrt{15}\alpha^{\frac{3}{2}}\tan^{-1}(\frac{8\sqrt{3}\pi
M}{\sqrt{5\alpha}})\Big{]},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{iii}^{TC}$ $\displaystyle=$
$\displaystyle-\frac{64\gamma\delta(4\pi)^{\delta+4}}{\sigma}\int_{M}^{0}\frac{m^{2+2\delta}dm}{(8\pi
m)^{2}+\frac{15}{9}\alpha},$ (19)

where hypergeometric functions are solutions to the first and third cases. It
is also useful to mention that, in all cases, the integrands recover Eq. (17)
at the corresponding appropriate limit. Since the third law is satisfied for
the first case and other cases when $0<\delta<\frac{1}{2}$, and $\delta<0$,
respectively, we plot the obtained evaporation times for some values of
$\delta$ that fall into these intervals.

Figure (2) shows the behavior of decay time against black hole mass for
different values of $\delta$ parameter. Each point on the curves presents the
time needed for a black hole of mass $0<M\leq 1$ to completely decay to zero
mass. As we observe within the upper panel, the decay time is finite and grows
as the initial black hole mass increases and asymptotically reaches a finite
value. Indeed, for $0<\delta<\frac{1}{2}$, it takes finite time for a black
hole to evaporate. The middle panel shows the behavior of decay time for
$\delta<1/2$ (family of black curves) and $1/2<\delta<1$ (family of red and
blue curves) where we observe that there exists a critical value of this
parameter ($\delta_{c}$) so that for $\delta<\delta_{c}$ the decay time is
finite and for $\delta>\delta_{c}$ the decay time diverges. In the lower
panel, we sketched the behavior of decay time for $\delta>1$ where it is seen
that the decay time for a black hole of finite mass grows unboundedly and
diverges as the initial mass of black hole increases. However, for these
values of $\delta$ parameter and those of middle panel, the third law of
thermodynamics ($S(T\rightarrow 0)\rightarrow 0$) is violated.

In Fig.(3) we have plotted decay time of a black hole against its mass for the
second case, where we observe that the larger the black hole mass the longer
it takes for the black hole to completely evaporate. The slope of
$\tau^{TC}_{ii}$ curve increases for larger values of black hole mass so that
the decay time grows unboundedly and diverges for massive and super massive
black holes. Finally, Fig. (4) presents the decay time of black hole for the
third case. In the upper panel we observe that, for $\delta<0$ for which the
third law is respected, a black hole with initial finite mass will completely
evaporate at a finite amount of time. The lower panel shows that for
$\delta>0$ the decay time is an increasing function of the black hole mass and
the heavier the initial black hole the longer it takes to completely
evaporate. However, from the viewpoint of the third case, for this value of
$\delta$ parameter the third law is not respected. As yet there is no
agreement on the numeric value of $\alpha$ parameter we have considered the
value of this parameter to be unity cqg ; agha . We further note that entropy
is a dimensionless quantity in each system of units where the Boltzmann
constant is unity. Hence, from Eq. (2) we can deduce that $\gamma\propto
C/\ell_{\rm Pl}^{2\delta}$ and as we work in the system of units for which
$\ell_{\rm Pl}=1$, $\gamma$ parameter is a positive constant value which we
have considered it to be unity.

Figure 2: Plot of decay time for the first case, versus black hole mass for
different values of parameter $\delta$. We have set $\sigma=\pi^{2}/60$,
$\gamma=1$ and $\alpha=1$. Figure 3: Plot of decay time for the second case,
versus black hole mass. We have set $\sigma=\pi^{2}/60$ and $\alpha=1$.

Figure 4: Plot of decay time for the third case, versus black hole mass. We
have set $\sigma=\pi^{2}/60$, $\gamma=1$ and $\alpha=1$.

###  V.2 Tsallis black hole

Although Eq. (16) is used in some previous works which investigate the black
hole thermodynamics in various non-extensive statistics epl ; refgerg ; gerg ;
mesri2 , a comprehensive study in the Tsallis framework should employs the
Tsallis counterpart of Eq. (16). The latter is a controversial issue tbbr1 ;
tbbr2 ; tbbr3 ; kol , as different averaging methods are useable and are
employed in this statistics kol . These different methods have their own
benefits and shortcomings, and indeed, their correctness and accessibility
situations are still unsolved and need more attentions and observations kol .
Here, motivated by the fact that

$\displaystyle\frac{dE}{dt}=-A\sigma_{q}T^{4}$ (20)

is obtained by different approaches tbbr1 , and is also originated from the
black body spectrum in Tsallis statistics tbbs , we focus on Eq. (20) as the
alternative of Eq. (16). Here, $\sigma_{q}$ is called the generalized SB
constant calculated as tbbr1

$\displaystyle\sigma_{q}=\frac{1}{\pi^{2}}\int_{0}^{\infty}\left[\frac{x^{3}}{\exp(x)-1}-\frac{1-q}{2}\frac{x^{5}\exp(x)}{(\exp(x)-1)^{2}}\right]dx\Rightarrow$
$\displaystyle\sigma_{q}\simeq\sigma(1-6\cdot 15(1-q))=\sigma(6\cdot
15q-5\cdot 15),$ (21)

if the integration is numerically solved tbbr1 . It is also obvious that
$\sigma_{q}\rightarrow\sigma$ at the appropriate limit of $q\rightarrow 1$.
For three obtained cases, we have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{i}^{q}=-\frac{(8\pi)^{3}}{2\sigma_{q}}\int_{M}^{0}m^{2}\exp\big{(}(1-q)4\pi
m^{2}\big{)}dm,$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{ii}^{q}=\frac{(8\pi)^{3}}{2\sigma_{q}}\frac{M^{3}}{3},$
(22) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{iii}^{q}=\tau_{i}^{q}.$

Since only the second and third cases satisfy the third law only for $q>1$, in
Fig. (5), we plotted $\tau_{ii}^{q}$ and $\tau_{iii}^{q}$ for different values
of $q$ parameter (family of black curves). We therefore observe that the decay
time is finite and the lesser the value of $q$ parameter, the longer it takes
for a Tsallis black hole to completely evaporate. The blue curve presents the
behavior of $\tau_{ii}^{q}$ where we see that the decay time grows unboundedly
as the initial black hole mass increases.

Figure 5: Plot of decay time for Tsallis black hole against its initial mass.
We have set $\sigma=\pi^{2}/60$. Figure 6: Plot of decay time for Kaniadakis
black hole against its initial mass. We have set $\sigma=\pi^{2}/60$.

###  V.3 Kaniadakis black hole

Black body spectrum in Kaniadakis statistics has recently been studied kanbb1
; kanbb2 , and it has been shown that kanbb1

$\displaystyle\frac{dE}{dt}=-A\sigma_{\kappa}T^{4},$ (23)

where $\sigma_{\kappa}=\frac{J_{3}^{\kappa}(0)}{4\pi^{2}}$ in which
$J_{3}^{\kappa}(0)=\int_{0}^{\infty}\frac{x^{3}}{\exp_{\kappa}(x)-1}dx$ while
$\exp_{\kappa}(x)=[\sqrt{1+\kappa^{2}x^{2}}+\kappa x]^{\frac{1}{\kappa}}$, and
we have $\sigma_{\kappa\rightarrow 0}=\sigma$ kanbb1 . Finally, we reach

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{i}^{\kappa}=-\frac{(8\pi)^{3}}{2\sigma_{\kappa}}\int_{M}^{0}m^{2}\cosh\big{(}4\kappa\pi
m^{2}\big{)}dm,$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{ii}^{\kappa}=\frac{(8\pi)^{3}}{2\sigma_{\kappa}}\frac{M^{3}}{3},$
(24)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\tau_{iii}^{\kappa}=\tau_{i}^{\kappa},$

for three cases we discussed above. In Fig. (6) we have plotted
$\tau_{ii}^{\kappa}$ and $\tau_{iii}^{\kappa}$ (family of black curves) for
some values of $\kappa<0$ parameter where we observe that the decay time grows
as $\kappa$ tends to larger values in negative direction. Such a behavior is
parallel to the satisfaction of the third law. For the second case we observe
that the decay time is finite for a finite mass black hole.

##  VI Summary

Based on the third law of thermodynamics, it is impossible for a system to
touch zero-entropy state (at least a state with minimum and finite value of
entropy) as its temperature tends to zero by only experiencing a finite number
of thermodynamical processes. As each process spends its own time interval to
be completed, we may even swap finite number of thermodynamical processes with
finite time in the above statement. The story becomes more complicated in the
ordinary black hole physics, generated by Bekenstein entropy, where it is
obtained that entropy diverges, while its temperature approaches zero. In this
situation, while black hole evaporates at finite time (17), and loses all of
its mass, its temperature diverges at its final evolution steps meaning that
we face a catastrophe ali . Here, we only focused on Schwarzschild black hole
as a primary solution which clearly exposes the mentioned inconsistency with
third law.

Motivated by the long range nature of gravity, some recent works that propose
a deep connection between quantum gravity and generalized statistics epl ;
homa ; mesri ; barrow ; mesri2 , and also the successes of these types of
statistics in justifying some cosmological and gravitational phenomena tsallis
; refgerg ; gerg ; non13 ; nonK ; EPJC ; KHDE ; sadeghi ; mesri2 , we studied
the status of third law for a Schwarzschild black hole in the framework of
some generalized statistics and we found out that this law may theoretically
be settled. Moreover, we obtained that the thermodynamic analysis along with
the laws may help us find new energy definitions, thus establishing a
consistency between the results of thermodynamics and the predictions of
quantum field theory about the black body radiation. The latter means that
thermodynamics may eventually shed light on the physics and pave the way to
find a comprehensive energy definition energy ; energy1 .

It is finally useful to mention that, a few months after submitting our
preprint to arXiv, we found two related papers plb ; gen1 . While one of them
investigates the availability of the generalized second law of thermodynamics
in a universe whose apparent horizon meets the generalized entropies plb ,
another one studies the BH temperature and energy by employing the Tsallis and
Cirto and Rényi entropies gen1 . The approach of Ref. gen1 , and additionally,
its findings about the Tsallis and Cirto entropy are similar to what we did
and obtained in Sec. (II), considered as a confirmation for our concern and
the strategy adopted here.

###### Acknowledgements.

IPL would like to acknowledge the contribution of the COST Action CA18108. IPL
was partially supported by the National Council for Scientific and
Technological Development - CNPq grant 306414/2020-1. J. P. M. G is supported
by CNPq under Grant No. 151701/2020-2.

## References

  * (1) I. Rácz, Class. Quantum Grav. 17, 4353 (2000).
  * (2) D. Chen et al, Int. J. Mod. Phys. A 29, 1430054 (2014).
  * (3) W. Xu, J. Wang, X. h. Meng, Galaxies 3, 53 (2015).
  * (4) Y. Yao, M. S. Hou, Y. C. Ong, Eur. Phys. J. C 79, 513 (2019).
  * (5) M. Hajebrahimi, K. Nozari, 4, 043E03 (2020).
  * (6) C. Tsallis, Introduction to Non-Extensive Statistical Mechanics: Approaching a Complex World, Springer, Berlin (2009).
  * (7) M. Masi, Phys. Lett. A 338, 217 (2005).
  * (8) C. Tsallis, L. J. L. Cirto, Eur. Phys. J. C 73, 2487 (2013).
  * (9) H. Moradpour, Int. Jour. Theor. Phys. 55, 4176 (2016).
  * (10) R. C. Nunes et al, JCAP 08, 051 (2016).
  * (11) N. Komatsu, Eur. Phys. J. C 77, 229 (2017).
  * (12) J. Sadeghi, M. Rostami, M. R. Alipour, Int. J. Mod. Phys. A 34 (30), 1950182 (2019).
  * (13) H. Moradpour et al, Eur. Phys. J. C 80, 732 (2020).
  * (14) T. S. Biró, V. G. Czinner, Phys. Lett. B 726, 861 (2013).
  * (15) S. Ghaffari et al, Gen. Rel. Grav. 51, 93 (2019).
  * (16) K. Mejrhit, R. Hajji, Eur. Phys. J. C 80, 1060 (2020).
  * (17) H. Moradpour et al, EPL 127, 60006 (2019).
  * (18) H. Shababi, K. Ourabah, Eur. Phys. J. Plus 135, 697 (2020).
  * (19) K. Mejrhit, S. E. Ennadifi, Phys. Lett. B 794, 24 (2019).
  * (20) J. D. Barrow, Phys. Lett. B 808, 135643 (2020).
  * (21) A. V. Kolesnichenko, Sol. Syst. Res. 54, 420 (2020).
  * (22) S. Nojiri, S. D. Odintsov, V. Faraoni, arXiv:2201.02424.
  * (23) S. W. Hawking, Nature. 248, 5443 (1974).
  * (24) J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * (25) A. Coutant, arXiv:1405.3466 [hep-th].
  * (26) H. Moradpour et al, AHEP, Article ID 7124730, (2018).
  * (27) M. Cruz, F. Izaurieta, S. Lepe, Eur. Phys. J. C 80, 559 (2020).
  * (28) L. C. Andrews, Special functions of mathematics for engineers, SPIE Press. p. 110 (1998).
  * (29) M. Cavaglia, Phys. Lett. B 569, 7 (2003).
  * (30) A. F. Ali, Phys. Rev. D 89, 104040 (2014).
  * (31) I. Lourek, M. Tribeche, Phys. Lett. A 381, 452 (2017).
  * (32) R. E. Kelly, Am. J. Phys. 49, 714 (1981).
  * (33) K. Nozari, S. F. Anvari, A. S. Sefiedgar, arXiv:1206.5631 [hep-th].
  * (34) V. Husain, S. S. Seahra, E. J. Webster, Phys. Rev. D 88, 024014 (2013).
  * (35) G. Amelino-Camelia, M. Arzano, Y. Ling, G. Mandanici, Class. Quant. Grav. 23, 2585 (2006).
  * (36) I. P. Lobo, G. B. Santos, Phys. Lett. B 817, 136272 (2021).
  * (37) S. Aghababaei et al, Phys. Scr. 96 055303 (2021).
  * (38) A. R. Plastino, A. Plastino, H. Vucetich, Phys. Lett. A 207, 42 (1995);
Q. A. Wang, A. Le Méhauté, Phys. Lett. A 237, 28 (1997);
U. Tirnakli, F. Büyükkiliç, D. Demirhan, Phys. Lett. A 245, 62 (1998).

  * (39) E. K. Lenzi, R. S. Mendes, Phys. Lett. A 250, 270 (1998).
  * (40) S. L. Choudhury, R. K. Paul, Annals of Physics, 395, 317 (2018).
  * (41) C. Tsallis, F. C. Sa Barreto, Generalization of Planck radiation law: possible testing of the nature of space-time, CBPF preprint. Brazil.
  * (42) K. Ourabah, M. Tribeche, Phys. Rev. E 89, 062130 (2014).
  * (43) S. Nojiri, S. D. Odintsov, V. Faraoni, Phys. Rev. D 104, 084030 (2021).
  * (44) E. M. C. Abreu, J. A. Neto, Phys. Lett. B 824, 136803 (2022).


# Veech groups, irrational billiards and stable abelian differentials

Ferrán Valdez  Max Planck Institut für Mathematik Vivatsgasse 7. 53111, Bonn,
Germany<EMAIL_ADDRESS>
###### Abstract.

We describe Veech groups of flat surfaces arising from irrational angled
polygonal billiards or irreducible stable abelian differentials. For
irrational polygonal billiards, we prove that these groups are non-discrete
subgroups of $\rm SO(2,\mathbf{R})$ and we calculate their rank.

##  1\. Introduction

The Veech group of a flat surface is the group of derivatives of orientation-
preserving affine homeomorphisms. If the surface is compact, Veech groups are
discrete subgroups of $\mathbf{SL}(2,\mathbf{R})$ that can be related to the
geodesic flow on the surface [7]. Our main goal is to describe Veech groups
arising from non-compact flat surfaces associated to billiards on an
irrational angled polygon. Nevertheless, in this article we will not discuss
dynamical aspects of geodesics. More precisely,

######  Theorem 1.1.

Let $P\subset\mathbf{R}^{2}$ be a simply connected polygon with interior
angles $\\{\lambda_{j}\pi\\}_{j=1}^{N}$, $S(P)$ the flat surface obtained from
$P$ via the Katok-Zemljakov construction and $G(S)$ the Veech group of $S(P)$.
Suppose there exists $\lambda_{j}\in\mathbf{R}\setminus\mathbf{Q}$ for some
$j=1,\ldots,n$. Then, $G(S)<\mathbf{SO}(2,\mathbf{R})$ and the group generated
by the rotations

(1.1) $R(S)=<\begin{pmatrix}\cos(2\lambda_{j}\pi)&-\sin(2\lambda_{j}\pi)\\\
\sin(2\lambda_{j}\pi)&\cos(2\lambda_{j}\pi)\end{pmatrix}\mid j=1,\ldots,N>$

has maximal rank in $G(S)$.

The surface $S(P)$ has infinite genus and only one end. A topological surface
satisfying these two conditions is called a _Loch Ness monster_ [6].
After a suggestion from M. Möller, we consider Veech groups arising from
stable abelian differentials at $\partial\Omega\overline{M_{g}}$, the boundary
of the Deligne-Mumford compactification of the Hodge bundle $\Omega M_{g}$
(see §5, [1] for a definition). On this boundary, the notion of Veech group
makes sense only for stable abelian differentials on an irreducible stable
curve, called irreducible. In this direction, we prove the following:

######  Proposition 1.2.

Let $(X,\omega)\in\partial\rm\Omega\overline{\mathcal{M}_{g}}$ be an
irreducible stable Abelian differential of genus $g$. Suppose that there
exists at least one node in $X$ where the 1–form $\omega$ has a pole. Let
$\\{r_{j},-r_{j}\\}_{j=1}^{k}$ be the set of all residues of $\omega$ and
define

(1.2) $\ N:=<\\{\begin{pmatrix}1&s\\\ 0&t\end{pmatrix}\mid
t\in\mathbf{R}^{+},s\in\mathbf{R}\\},-Id>.$

Let $G(X)=G(X,\omega)$ be the Veech group of $(X,\omega)$. Then,

  1. (1)

If there exist $i\neq j$ such that $r_{i}/r_{j}\notin\mathbf{R}$, then $G(X)$
is finite.

  2. (2)

If all residues of $\omega$, as vectors in $\mathbf{C}\simeq\mathbf{R}^{2}$
are parallel, then $G(X)<N$ is either conjugated to a discrete subgroup or the
equal to $N$.

Recently, Hubert and Schmithüsen [3] have shown the existence of a countably
family of _infinite area_ origamis whose Veech groups are infinitely generated
subgroups of $\mathbf{SL}(2,\mathbf{Z})$. These origamis arise as
$\mathbf{Z}$-covers of (finite area) genus 2 origamis. Motivated by this work,
in the last section of this article we construct, for each $n\in\mathbf{N}$,
an uncountable family of flat surfaces $\mathcal{S}_{n}=\\{S_{i}\\}_{i\in I}$
such that each $S_{i}$ is homeomorphic to the Loch Ness monster and the Veech
group $G(S_{i})<\mathbf{SO}(2,\mathbf{R})$ is infinitely generated.

This article is organized as follows. We introduce the notion of _tame_ flat
surface and extend the definition of some classical geometric invariants
(saddle connections, Veech groups) to the non-compact realm in Section 2 .
Loosely speaking, _tame_ flat surfaces present a discrete set of
singularities, which are either of finite or infinite angle. We briefly recall
the Katok-Zemljakov construction, the notion of stable Abelian differential
and define Veech groups for irreducible nodal flat surfaces. Section 3 deals
with the proof of Theorem 1.1 and Section 4 with the proof of Proposition 1.2.
Finally, Section 5 presents the construction of the family of flat surfaces
$\mathcal{S}_{n}$ mentioned above.

Acknowledgments. This article was written during a stay of the author at the
Max Planck Institut für Mathematik in Bonn. The author wishes to express his
gratitude to the administration and staff of the MPI for the wonderful working
facilities and the atmosphere. The author acknowledges support from the
Sonderforschungsbereich/Transregio 45 and the ANR Symplexe. The author thanks
M. Möller and M. Bainbridge for valuable discussions.

##  2\. Preliminaries

Non-compact flat surfaces. Let $(S,\omega)$ be a pair formed by a connected
Riemann surface $S$ and a holomorphic 1–form $\omega$ on $S$ which is not
identically zero. Denote by $Z(\omega)\subset S$ the zero locus of the form
$\omega$. Local integration of this form endows $S\setminus Z(\omega)$ with an
atlas whose transition functions are translations of $\mathbf{C}$. The
pullback of the standard translation invariant flat metric on the complex
plane defines a flat metric $d$ on $S\setminus Z(\omega)$. Let $\widehat{S}$
be the metric completion of $S$. Each point in $Z(\omega)$ has a neighborhood
isometric to the neighborhood of $0\in\mathbf{C}$ with the metric induced by
the 1–form $z^{k}dz$ for some $k>1$ (which is a cyclic finite branched
covering of $\mathbf{C}$). The points in $Z(\omega)$ are called finite angle
singularities. Note that there is a natural embedding of $S$ into
$\widehat{S}$.

######  Definition 2.1.

A point $p\in\widehat{S}$ is called an _infinite angle singularity_ , if there
exists a radius $\epsilon>0$ such that the punctured neighborhood:

(2.3) $\\{z\in\widehat{S}\mid 0<d(z,p)<\epsilon\\}$

is isometric to the infinite cyclic covering of
$\epsilon\mathbf{D}^{*}=\\{w\in\mathbf{C}^{*}\mid 0<\mid w\mid<\epsilon\\}$.
We denote by $Y_{\infty}(\omega)$ the set of infinite angle singularities of
$\widehat{S}$.

######  Definition 2.2.

The pair $(S,\omega)$ is called a _tame_ flat surface if $\hat{S}\setminus
S=Y_{\infty}(\omega)$.

One can easily check that flat surfaces arising from irrational polygons or
stable abelian differentials are tame.

######  Definition 2.3.

A _singular geodesic_ of $S=(S,\omega)$ is an open geodesic segment in the
flat metric $d$ whose image under the natural embedding
$S\hookrightarrow\widehat{S}$ issues from a singularity of $\widehat{S}$,
contains no point of $Y(\omega)$ in its interior and is not properly contained
in some other geodesic segment. A _saddle connection_ is a finite length
singular geodesic.

To each saddle connection we can associate a _holonomy vector_ : we ’develop’
the saddle connection in the plane by using local coordinates of the flat
structure. The difference vector defined by the planar line segment is the
holonomy vector. Two saddle connections are _parallel_ , if their
corresponding holonomy vectors are linearly dependent.
Let $\mathrm{Aff}_{+}(S)$ be the group of affine orientation preserving
homeomorphisms of the flat surface $S$ (by definition $S$ comes with a
distinguished 1–form $\omega$). Consider the map

(2.4)
$\mathrm{Aff}_{+}(S)\overset{D}{\longrightarrow}\mathbf{GL}_{+}(2,\mathbf{R})$

that associates to every $\phi\in\mathrm{Aff}_{+}(S)$ its (constant) Jacobian
derivative $D\phi$.

######  Definition 2.4.

Let $S$ be a flat surface. We call $G(S)=D(\mathrm{Aff}_{+}(S))$ the _Veech
group_ of $S$.

The Katok-Zemljakov construction. In the following paragraph we recall briefly
the definition of this construction. For details see [6] and references
within.
Let $P_{0}$ denote the polygon $P$ deprived of its vertices. The
identification of two disjoint copies of $P_{0}$ along ”common sides” defines
a Euclidean structure on the $N$-punctured sphere. We denote it by
$\SS^{2}(P)$. This punctured surface is naturally covered by $S(P_{0})$, the
_minimal translation surface_ corresponding to $P$. We denote the projection
of this covering by $\pi:S(P_{0})\longrightarrow\SS^{2}(P)$. Call a vertex of
$P$ _rational_ , if the corresponding interior angle is commensurable with
$\pi$. When the set of rational vertices of $P$ is not empty, the translation
surface $S(P_{0})$ can be locally compactified by adding points ”above”
rational vertices of $P$. The result of this local compactification is a flat
surface with finite angle singularities that we denote by $S(P)$. If the set
of rational vertices of $P$ is empty, we set $S(P)=S(P_{0})$. In both cases,
$S(P)$ is called the flat surface obtained from the polygon $P$ via the
_Katok-Zemljakov construction_. Remark that, in the case of rational polygons,
some authors give a different definition (see [5] or [2]).

Stable Abelian differentials. We recall briefly the notion of stable Abelian
differential, following Bainbridge [1].
A _nodal Riemann surface_ $X$ is a finite type Riemann surface, _i.e._ with
finitely generated fundamental group, that has finitely many cusps which have
been identified pairwise to form nodes. A connected component of a nodal
Riemann surface $X$ with its nodes removed is called a _part_ of $X$, and the
closure of a part of $X$ is an _irreducible component_. The genus of a nodal
Riemann surface is the topological genus of the non-singular Riemann surface
obtained by replacing each node in $X$ with an annulus. A _stable Riemann
surface_ is a connected nodal Riemann surface for which each part has negative
Euler characteristic. A _stable Abeliann differential_ $\omega$ on a stable
Riemann surface $X$ is a holomophic 1–form on $X$ minus its nodes such that
its restriction to each part of $X$ has at worst simple poles at the cusps and
at two cusps which have been identified to form a node, the differential has
opposite residues, if any. Nodes at which $\omega$ presents a pole are called
_polar nodes_.

Veech groups on stable Abelian differentials. Let $(X,\omega)$ be a stable
Abelian differential. We denote by ${\rm Aff_{+}}(X,\omega)$ the group of
affine orientation preserving homeomorphisms of $X$. The Jacobian derivative
$D\phi$ of an affine homeomorphism $\phi$ is constant on each irreducible
component of $(X,\omega)$. In general, there is no canonical derivation
morphism from the affine group of a stable Abelian differential onto
$\mathbf{GL}_{+}(2,\mathbf{R})$. Consider, for example, the genus 2 stable
Abelian differential given by the following figure:

Figure 1.

We avoid this situation by restricting ourselves to irreducible Riemann
surfaces.

######  Definition 2.5.

Let $X=(X,\omega)$ be an _irreducible_ stable Abelian differential. We call
$G(X)=D({\rm Aff_{+}}(X,\omega))$ the _Veech group_ of $X$.

Abelian differentials close to a stable Abelian differential $\rm(X,\omega)$
with a polar node develop very long cylinders which are pinched off to form a
node in the limit, (see §5.3 [1]). In the following figure we depict a genus
two stable abelian differential with two nodes ( with residues $\rm\pm 1$ and
$\rm\pm(1+i)$) and two double zeroes:

Figure 2.

When considering the flat metric, every stable Abelian differential deprived
of its polar nodes is a complete metric space. We call _singular geodesic_ in
the context of stable Abelian differentials, every geodesic segment that
issues from a zero or a non-polar node of $\omega$, contains no such zero or
non-polar node on its interior and is not properly contained in some other
geodesic segment. As before, finite length singular geodesics will be called
_saddle connections_.

_Decomposition of stable Abelian differentials with polar nodes_. Suppose that
$(X,\omega)$ has polar nodes with residues $r_{1},\ldots,r_{k}$. Every $r_{j}$
defines a direction $\theta(r_{j})\in\mathbf{R}/\mathbf{Z}$ for which
$(X,\omega)$ presents a set of disjoint infinite area cylinders
$C_{1,j},\ldots,C_{n(j),j}$ foliated by closed geodesics parallel to
$\theta(r_{j})$ and whose length is $\mid r_{j}\mid$. Denote by $C_{j}$ the
closure in $(X,\omega)$ of $\cup_{i=1}^{n(j)}C_{i,j}$ and
$C=\cup_{j=1}^{k}C_{j}$. We define

(2.5) $X^{\prime}:=X\setminus C$

The Veech group of $(X,\omega)$ acts linearly on the set of residues of
$\omega$ and leaves the decomposition $X=X^{\prime}\sqcup C$ invariant.

##  3\. Proof of Theorem 1.1

First, we prove that the matrix group $R(S)$ defined in (1.1) is a subgroup of
$G(S)$. Then, we prove that $G(S)<\mathbf{SO}(2,\mathbf{R})$ and, finally,
that ${\rm Rank}(G(S))={\rm Rank}(R(S))$.

(i) The locally Euclidean structure on the $N$-punctured sphere
$\mathbf{S}^{2}(P)$ gives rise to the holonomy representation:

(3.6) $\rm hol:\pi_{1}(\mathbf{S}^{2}(P))\longrightarrow
Isom_{+}(\mathbf{R}^{2})$

Let $B_{j}$ be a simple loop in $\mathbf{S}^{2}(P)$ around the missing vertex
of $P$ whose interior angle is $\lambda_{j}\pi$, $\rm j=1,\ldots,N$. Suppose
that $B_{j}\cap B_{i}=*$, for $\rm i\neq j$. Then, $\\{B_{j}\\}_{j=1}^{N}$
generates $\pi_{1}(\mathbf{S}^{2}(P),*)$. Given an isometry $\varphi\in\rm
Isom_{+}(\mathbf{R}^{2})$, we denote its derivative by $D\circ\varphi$. A
direct calculation in local coordinates shows that $\rm hol(B_{j})$ is affine
and that $M_{j}=D\circ\rm hol(B_{j})$ is given by:

(3.7) $M_{j}=\begin{pmatrix}\cos(2\lambda_{j}\pi)&-\sin(2\lambda_{j}\pi)\\\
\sin(2\lambda_{j}\pi)&\cos(2\lambda_{j}\pi)\end{pmatrix}\hskip
28.45274ptj=1,\ldots,N.$

Since $G(S(P_{0}))=G(S(P))$, we conclude that $R(S)$ is a subgroup of $G(S)$.
(ii) We claim that length of every saddle connection in $S(P)$ is bounded
below by some constant $c=c(P)>0$. Indeed, consider the folding map
$f:\SS^{2}(P)\longrightarrow P$ which is 2-1 except along the boundary of $P$.
The projection $f\circ\pi:S(P_{0})\longrightarrow P$ maps every saddle
connection $\gamma\subset S(P_{0})$ onto a _generalized diagonal_ of the
billiard game on $P$ (see [4] for a precise definition). The length of
$\gamma$ is bounded below by the length of the generalized diagonal
$f\circ\pi(\gamma)$. The length of any generalized diagonal of the billiard
table $P$ is bounded below by some positive constant $c$ depending only on
$P$. This proves our claim. The constant $c$ is realized by a generalized
diagonal. Therefore, we can choose a holonomy vector $v$ is of minimal length.
Given that $R(S)<G(S)$, the $G(S)$-orbit of $v$ is dense in the circle of
radius $|v|$ centered at the origin. This forces the Veech group $G(S)$ to lie
in $\mathbf{SO}(2,\mathbf{R})$.

(iii) Suppose that there exist an affine homeomorphism $\varphi\in{\rm
Aff}_{+}(S)$ such that $D\varphi$ is an infinite order element of
$\mathbf{SO}(2,\mathbf{R})/R(S)$. Let $\gamma_{0}$ be a fixed saddle
connection. Then $\\{f\circ\pi\circ\varphi^{k}(\gamma_{0})\\}_{k\in Z}$ is an
infinite set of generalized diagonals of bounded length. But this is a
contradiction, for the set of generalized diagonals of bounded length on a
polygonal billiard is always finite [4].

##  4\. Proof of Proposition 1.2

The Veech group of the irrational stable Abelian differential $(X,\omega)$
acts linearly on the (finite) set of residues of $\omega$. Therefore, if not
all residues are parallel, $G(X)$ must be finite.
Suppose now that all residues are parallel to the horizontal direction. Then
$G(X)<N$. If every holonomy vector of $(X,\omega)$ is horizontal we claim that
$G(X)=N$. Indeed, in this situation $X^{\prime}$ defined in (2.5) is empty and
the horizontal geodesic flow decomposes $X$ into finitely many cylinders with
horizontal boundaries. This allows to define, for every $g\in N$, an
orientation preserving affine homeomorphism of $X$ whose differential is
exactly $g$. On the other hand, if at least one holonomy vector fails to be
horizontal, then $G(X)<N$ is discrete, for the set of holonomy vectors of any
stable Abelian differential is discrete.
$\square$

Remark. Veech groups of irreducible stable Abelian differentials in
$\partial\Omega\overline{\mathcal{M}_{g}}$ without polar nodes are as
”complicated” as Veech groups of flat surfaces in $\Omega\mathcal{M}_{g}$ with
marked points. More precisely, a nodal Riemann surface $X$ has a
_normalization_ $f:S(X)\longrightarrow X$ defined by separating the two
branches passing through each node of $X$. For every node $p$, denote
$\\{p_{+},p_{-}\\}:=f^{-1}(p)$. Then, if the stable Abelian differential
$(X,\omega)$ has no polar nodes, we have the equality:

(4.8) ${\rm Aff}_{+}(X,\omega)=\\{\phi\in{\rm Aff}_{+}(S(X),\omega)\hskip
2.84526pt|\hskip 2.84526pt\phi(p_{+})=\phi(p_{-}),\hskip
2.84526pt\forall\hskip 2.84526ptp\hskip 2.84526pt\text{node of $X$}\\}.$

##  5\. Infinitely generated Veech groups in $\mathbf{SO}(2,\mathbf{R})$

Fix $n\in\mathbf{N}$. Consider an unbounded sequence of real numbers

(5.9) $\rm x_{0}=0<x_{1}<x_{2}<\ldots<x_{j}<...$

such that $\rm x_{j+1}-x_{j}>1$ for all $j$. The segments of straight line
joining the point $\rm(x_{j},x_{j}^{2n})$ to$\rm(x_{j+1},x_{j+1}^{2n})$ and
$\rm(-x_{j},x_{j}^{2n})$ to $\rm(-x_{j+1},x_{j+1}^{2n})$, $j\geq 0$, define a
polygonal line $\partial P$ in $\mathbf{C}$. Let $int(P)$ be the connected
component of $\mathbf{C}\setminus\partial P$ intersecting the positive
imaginary axis $\rm Im(z)>0$. We define $P=\partial P\cup int(P)$. We call $P$
the _unbouded polygon_ defined by the sequence (5.9). Remark that $P$ is
symmetric with respect to the imaginary axis. For each $\rm j\geq 0$, let
$\rm\lambda_{j}\pi$ be the interior angle of $P$ at the vertex
$\rm(x_{j},x_{j}^{2n})$.

######  Definition 5.1.

We say that a sequence of real numbers $\rm\\{\mu_{j}\\}_{j\geq 0}$ is _free
of resonances_ if and only if for every finite subset
$\rm\\{\mu_{j_{1}},\ldots\mu_{j_{N}}\\}$ the kernel of the group morphism
$\mathbf{Z}^{N}\longrightarrow\mathbf{C}$ defined by

$\rm(n_{1},\ldots,n_{N})\longrightarrow exp(2\pi
i(\sum_{k=1}^{N}n_{k}\mu_{j_{k}}))$

is trivial.

There are uncountable many choices $\\{x_{j\geq 0}^{i}\\}$, $i\in I$, for
(5.9), such that the sequence $\\{\lambda_{j}^{i}\\}_{j\geq 0}$ defining the
interior angles of $P=P_{i}$ is free of resonances. For each $i\in I$, denote
by $S^{2}(P_{i})$ the identification of two vertexless copies of $P_{i}$ along
”common sides”. The Katok-Zemljakov construction described in Section 2 can be
applied to the unbounded polygon $P_{i}$. The result is a flat covering
$S_{i}\longrightarrow\mathbf{S}^{2}(P_{i})$.

######  Lemma 5.2.

The flat surface $S_{i}$ is homeomorphic to the Loch Ness monster. The Veech
group $G(S_{i})<\mathbf{SO}(2,\mathbf{R})$ is infinitely generated and
contains the infinite rank group generated by the matrices

(5.10) $\begin{pmatrix}\cos(2\lambda_{j}^{i}\pi)&-\sin(2\lambda_{j}^{i}\pi)\\\
\sin(2\lambda_{j}^{i}\pi)&\cos(2\lambda_{j}^{i}\pi)\end{pmatrix}\rm\hskip
28.45274ptj\geq 0,$

###### Proof.

Every flat surface obtained via the Katok-Zemljakov construction from a
(bounded) polygon whose angles are all irrational multiples of $\pi$, is
homeomorphic to a Loch-Ness monster. This is proved in Theorem 1 ( Case (A)
and absence of resonances) in [6]. The same conclusion can be drawn for the
unbounded polygons $P_{i}$ after replacing in the proof of Theorem 1 [_Ibid._]
_polygon P_ and surface $X(P)$ by _unbounded polygon_ $P_{i}$ and $S_{i}$,
respectively.
In §3, sections (i) and (ii) made use of the boundness of $P$ to assure that
the length of every saddle connection in $S(P)$ was bounded below by a
constant depending only on $P$. For unbounded polygons, this is ensured by the
condition $\rm x_{j+1}-x_{j}>1$, for all $j$, on the sequence (5.9). It
follows that, for every $i\in I$, $G(S_{i})<\mathbf{SO}(2,\mathbf{R})$ and
that this Veech group contains the group generated by the matrices (5.10).
This matrix group is infinitely generated, for the sequence
$\\{\lambda_{j}^{i}\\}_{j\geq 0}$ is free of resonances, for every $i\in I$.
$\square$

## References

  * [1] M. Bainbridge. _Euler characteristics of Teichmüller curves in genus two_. Geom. Topol. 11 (2007), 1887–2073.
  * [2] E. Gutkin and S. Troubetzkoy _Directional flows and strong recurrence for polygonal billiards_. International Conference on Dynamical Systems (Montevideo, 1995), 21–45, Pitman Res. Notes Math. Ser., 362, Longman, Harlow, 1996. _Geom. Dedicata_ 125 (2007), 39–46.
  * [3] P. Hubert and G. Schmithüsen. _Infinite translation surfaces with infinitely generated Veech groups_. Preprint. http://www.cmi.univ-mrs.fr/ hubert/articles/hub-schmithuesen.pdf
  * [4] A. B. Katok, _The growth rate for the number of singular and periodic orbits for a polygonal billiard_. Comm. Math. Phys. 111 (1987), no. 1, 151–160.
  * [5] H. Masur and S. Tabachnikov. _Rational billiards and flat structures_. Handbook of dynamical systems. Vol. 1A, 1015–1089, North Holland. Amsterdam 2002.
  * [6] J.F. Valdez. _Infinite genus surfaces and irrational polygonal billiards_. To appear in Geom. Dedicata.
  * [7] W.A. Veech. _Teichmüller curves in the moduli space, Eisenstein series and applications to triangular billiards_. Inventiones mathematicae, 97, 1989, 553–583.