[table]capposition=top capbtabboxtable[][]

# Primordial non-Guassianity in inflation with gravitationally enhanced
friction

Li-Yang<EMAIL_ADDRESS>Hongwei Yu222<EMAIL_ADDRESS>and Puxun
Wu333<EMAIL_ADDRESS>Department of Physics and Synergetic Innovation Center
for Quantum Effects and Applications, Hunan Normal University, Changsha, Hunan
410081, China

###### Abstract

The gravitationally enhanced friction can reduce the speed of the inflaton to
realize an ultra-slow-roll inflation, which will amplify the curvature
perturbations. The amplified perturbations can generate a sizable amount of
primordial black holes (PBHs) and induce simultaneously a significant
background gravitational waves (SIGWs). In this paper, we investigate the
primordial non-Gaussianity of the curvature perturbations in the inflation
with gravitationally enhanced friction. We find that when the gravitationally
enhanced friction plays a role in the inflationary dynamics, the non-
Gaussianity is noticeably larger than that from the standard slow-roll
inflation. During the regime in which the power spectrum of the curvature
perturbations is around its peak, the non-Gaussianity parameter changes from
negative to positive. When the power spectrum is at its maximum, the non-
Gaussianity parameter is near zero ($\sim\mathcal{O}(0.01)$). Furthermore, the
primordial non-Gaussianity promotes the formation of PBHs, while its effect on
SIGWs is negligible.

##  I Introduction

Inflation resolves most of the problems, such as the flatness, horizon and
monopole problems, that plague the standard cosmological model Guth1980 ;
Linde1982 ; Starobinsky1980 ; Albrecht1982 . During inflation the curvature
perturbations are stretched outside the Hubble horizon and then stop
propagating with the amplitudes frozen at certain nonzero values. Inflation
predicts a nearly scale-invariant spectrum for the curvature perturbations,
which is well consistent with the CMB observations Aghanim2020 . The CMB
observations indicate that the amplitude $\mathcal{P_{R}}$ of the power
spectrum of the curvature perturbations is about $10^{-9}$ Aghanim2020 . After
inflation, these super-horizon perturbations, which will reenter the Hubble
radius during the radiation- or matter-dominated era, result in the formation
of large scale cosmic structures and at the same time lead to possible
generation of primordial black holes (PBHs) Hawking1971 ; Carr1974 ; Carr1975
; Khlopov2007 . The possibility is however slim for the standard slow-roll
inflation since the amplitude of the power spectrum of the curvature
perturbations is too small ($\sim 10^{-9}$).

If a sizable amount of PBHs is formed in the early universe, PBHs with
different masses can be used to explain different astronomical events. For
example, the $\mathcal{O}(10)M_{\odot}$, $\mathcal{O}(10^{-5})M_{\odot}$ and
$\mathcal{O}(10^{-12})M_{\odot}$ PBHs can explain the gravitational wave
events observed by the LIGO/Virgo collaboration lg1 ; lg2 ; lg3 ; lg4  and six
ultrashort-timescale microlensing events in the OGLE data P.Mroz2017 ;
H.Niikura2019 , and make up all dark matter A.Katz2018 ; H.Niikura20191 ;
A.Barnacka2012 ; P.W.Graham2015 ; Belotsky2014 , respectively, where
$M_{\odot}$ is the mass of the Sun. To generate abundant PBHs,
$\mathcal{P_{R}}$ is required to reach the order of $\mathcal{O}(10^{-2})$.
Since the CMB observations have put stringent constraints on $\mathcal{P_{R}}$
only at the CMB scales, we can realize the production of abundant PBHs by
enhancing the amplitude of the power spectrum of the curvature perturbations
about seven orders at small scales. As $\mathcal{P_{R}}\propto 1/\epsilon$
with $\epsilon$ being the slow roll parameter, a natural way to amplify the
curvature perturbations is to include an ultra-slow-roll period during
inflation. Flattening the inflationary potential can reduce the rolling speed
of the inflaton, which gives rises to an ultra-slow-roll inflation Germani2017
; Motohashi2017 ; Ezquiaga2017 ; H. Di2018 ; Ballesteros2018 ; Dalianis2019 ;
Gao2018 ; Drees2021 ; C.Fu2020 ; Xu2020 ; Lin2020 ; Dalianis2021 ; Yi2021 ;
Gao2021 ; Yi2021b ; TGao2021 ; Solbi2021 ; Gao2021b ; Solbi2021b ; Zheng2021 ;
Teimoori2021a ; Cai2021 ; Wang2021 ; Karam2022 . The ultra-slow-roll inflation
can also be achieved via slowing down the inflaton by gravitationally
enhancing friction Fuchengjie2019 ; fuchengjie2020 ; Dalianis(2020) ;
Teimoori2021 ; Heydari2022 ; Heydari2022b . Moreover, some other mechanisms,
such as parametric resonance yfcai2018 ; yfcai2019 ; c.chen2019 ; c.chen2020 ;
Addazi2022 ; Cai2020 , have also been proposed to amplify the curvature
perturbations.

When the amplified curvature perturbations reenter the Hubble horizon during
the radiation- or matter-dominated era, they will not only generate the PBHs,
but also lead simultaneously to large scalar metric perturbations, which
become an effective source of background gravitational waves. These
gravitational waves, called the scalar induced gravitational waves (SIGWs),
may be detectable by the future GW projects such as LISA lisa , Taiji taiji ,
TianQin tianqin  and PTA pta1 ; pta2 ; pta3 ; pta4 .

When we assess the abundance of PBHs and the energy density of SIGWs, the
curvature perturbations are assumed usually to be of a Gaussian distribution.
This is because the curvature perturbations generated during the standard
slow-roll inflation are nearly Gaussian with negligible non-Gaussianity.
However, once the inflation departs from the slow-roll inflation or it is
driven by the noncanonical fields, the primordial non-Gaussianity of the
curvature perturbations may no longer be ignored. The primordial non-
Gaussianity in the ultra-slow-roll inflation has been studied widely Cai2019 ;
QingGuoHuang2013 ; Zhang2021 ; F.Zhang2022 ; fengge2020 ; Chul-Moon Yoo2019 ;
G.Franciolini2018 ; Matthew2022 ; SamuelPassaglia2019 ; BravoRafael2018 ;
Cai2018 ; VicenteAtal2018 ; VicenteAtal2019 , becuase the abundance of PBHs is
extremely sensitive to the primordial non-Gaussianity of the curvature
perturbations. For the PBHs generated from inflation with gravitationally
enhanced friction mechanism Fuchengjie2019 ; Germani2011_1 ; Germani2011_2 ,
the primordial non-Gaussianity might be non-negligible too since the inflation
field couples derivatively with the gravity and the rolling of the inflaton is
ultra slow. In this paper we study, in the ultra-slow-roll inflation achieved
through gravitationally enhanced friction, the non-Gaussianity of the
curvature perturbations and its effect on the PBH abundance and the energy
density of SIGWs.

The paper is organized as follows: In Sec. II, we briefly review the inflation
model with the nonminimal derivative coupling between inflation field and
gravity. Sec. III studies the primordial non-Gaussianity of the curvature
perturbations. In Sec. IV, the effect of the non-Gaussianity of the curvature
perturbations on the abundance of PBHs and the energy density of SIGWs are
assessed. Finally, we give our conclusions in Sec. V.

##  II inflation with the gravitationally enhanced friction

To enhance the friction term in the equation of motion of the inflaton through
the gravity, we consider a nonminimal derivative coupling between the inflaton
field $\phi$ and gravity, with the action given by

$\displaystyle\mathcal{S}=\int
d^{4}x\sqrt{-g}\left[\frac{M_{\mathrm{pl}}^{2}}{2}R-\frac{1}{2}\left(g^{\mu\nu}-\frac{1}{M_{\mathrm{pl}}^{2}}\theta(\phi)G^{\mu\nu}\right)\nabla_{\mu}\phi\nabla_{\nu}\phi-V(\phi)\right],$
(1)

where $M_{\mathrm{pl}}$ is the reduced Planck mass, and $g$ is the determinant
of the metric tensor $g_{\mu\nu}$, $R$ is the Ricci scalar, $G_{\mu\nu}$ is
the Einstein tensor, $\theta(\phi)$ is the coupling function, and $V(\phi)$ is
the potential of the scalar inflaton field.

In the spatially flat Friedmann-Robertson-Walker background

$\displaystyle ds^{2}=-dt^{2}+a(t)^{2}d\mathbf{x}^{2}$ (2)

with $a(t)$ being the scale factor, one can obtain, from the action (1), the
background equations

$\displaystyle
3H^{2}=\frac{1}{M_{\mathrm{pl}}^{2}}\left[\frac{1}{2}\left(1+\frac{9}{M_{\mathrm{pl}}^{2}}\theta(\phi)H^{2}\right)\dot{\phi}^{2}+V(\phi)\right],$
(3)
$\displaystyle-2\dot{H}=\frac{1}{M_{\mathrm{pl}}^{2}}\left[\left(1+\frac{3}{M_{\mathrm{pl}}^{2}}\theta(\phi)H^{2}-\frac{1}{M_{\mathrm{pl}}^{2}}\theta(\phi)\dot{H}\right)\dot{\phi}^{2}-\frac{1}{M_{\mathrm{pl}}^{2}}\theta_{,\phi}H\dot{\phi}^{3}-\frac{2}{M_{\mathrm{pl}}^{2}}\theta(\phi)H\dot{\phi}\ddot{\phi}\right],$
(4)
$\displaystyle\left(1+\frac{3}{M_{\mathrm{pl}}^{2}}\theta(\phi)H^{2}\right)\ddot{\phi}+\left[1+\frac{1}{M_{\mathrm{pl}}^{2}}\theta(\phi)\left(2\dot{H}+3H^{2}\right)\right]3H\dot{\phi}+\frac{3}{2M_{\mathrm{pl}}^{2}}\theta_{,\phi}H^{2}\dot{\phi}^{2}+V_{,\phi}=0\,,$
(5)

where an overdot denotes the derivative with respective to the cosmic time
$t$, $H=\frac{\dot{a}}{a}$ is the Hubble parameter,
$\theta_{,\phi}=d\theta/d\phi$, and $V_{,\phi}=dV/d\phi$.

To describe the slow-roll inflation, we define the slow-roll parameters

$\displaystyle\epsilon$
$\displaystyle=-\frac{\dot{H}}{H^{2}},\quad\delta_{\phi}=\frac{\ddot{\phi}}{H\dot{\phi}},$
$\displaystyle\delta_{X}$
$\displaystyle=\frac{\dot{\phi}^{2}}{2{M_{\mathrm{pl}}^{2}}H^{2}},\quad\delta_{D}=\frac{\theta\dot{\phi}^{2}}{4{M_{\mathrm{pl}}^{4}}}.$
(6)

When $\\{\epsilon,|\delta_{\phi}|,\delta_{X},\delta_{D}\\}\ll 1$ are
satisfied, the slow-roll inflation is obtained.

Table 1: Three different parameter sets for generating the $\mathcal{O}(10)M_{\odot}$, $\mathcal{O}(10^{-5})M_{\odot}$ and $\mathcal{O}(10^{-12})M_{\odot}$ PBHs, respectively. $\phi_{*}$ and $N_{*}$ are the value of inflation field and e-folding number when the pivot scale $k_{*}=0.05\mathrm{Mpc}^{-1}$ exits the Hubble horizon.  | $\phi_{*}/M_{\mathrm{pl}}$ | $\phi_{c}/M_{\mathrm{pl}}$ | $\sigma_{s}$ | $\lambda$ | $w$ | m | $n_{s}$ | $r$ | $N_{*}$ |
---|---|---|---|---|---|---|---|---|---|---
Case 1 | $4.29$ | $4.02$ | $1.8\times 10^{-9}$ | $6.68\times 10^{-10}$ | $3.70\times 10^{16}$ | $6\times 10^{8}$ | 0.971 | 0.0350 | 64 |
Case 2 | $4.30$ | $3.63$ | $1.8\times 10^{-9}$ | $6.60\times 10^{-10}$ | $4.08\times 10^{16}$ | $8\times 10^{8}$ | 0.972 | 0.0340 | 66 |
Case 3 | $3.95$ | $2.95$ | $2.0\times 10^{-9}$ | $7.40\times 10^{-10}$ | $5.19\times 10^{16}$ | $9.5\times 10^{8}$ | 0.971 | 0.0357 | 68.5 |

In order to find the power spectrum of the curvature perturbations, we need to
derive the quadratic action for the curvature perturbations $\mathcal{R}$ from
the action given in Eq. (1), which takes the form A.D.Felice2011 ;
S.Tsujikawa2012 ; Kobayashi2011

$\displaystyle S^{(2)}=\int
dtd^{3}xa^{3}Q\left[\dot{\mathcal{R}}^{2}-\frac{c_{s}^{2}}{a^{2}}(\partial\mathcal{R})^{2}\right],$
(7)

where

$\displaystyle Q=\frac{w_{1}\left(4w_{1}w_{3}+9w_{2}^{2}\right)}{3w_{2}^{2}},$
(8) $\displaystyle
c_{s}^{2}=\frac{3\left(2w_{1}^{2}w_{2}H-w_{2}^{2}w_{4}+4w_{1}\dot{w}_{1}w_{2}-2w_{1}^{2}\dot{w}_{2}\right)}{w_{1}\left(4w_{1}w_{3}+9w_{2}^{2}\right)},$
(9)

and

$\displaystyle w_{1}$
$\displaystyle=M_{\mathrm{pl}}^{2}\left(1-2\delta_{D}\right),$ $\displaystyle
w_{2}$ $\displaystyle=2HM_{\mathrm{pl}}^{2}\left(1-6\delta_{D}\right),$
$\displaystyle w_{3}$
$\displaystyle=-3H^{2}M_{\mathrm{pl}}^{2}\left(3-\delta_{X}-36\delta_{D}\right),$
$\displaystyle w_{4}$
$\displaystyle=M_{\mathrm{pl}}^{2}\left(1+2\delta_{D}\right).$ (10)

From Eq. (7), we obtain the Mukhanov-Sasaki equation

$\displaystyle
u_{k}^{\prime\prime}+\left(c_{s}^{2}k^{2}-\frac{z^{\prime\prime}}{z}\right)u_{k}=0,$
(11)

where $z^{2}=2a^{2}Q$, and $u_{k}=z\mathcal{R}_{k}$. Solving this Mukhanov-
Sasaki equation yields the power spectrum of the curvature perturbations

$\displaystyle\mathcal{P}_{\mathcal{R}}\simeq\mathcal{P}_{\mathcal{R}{{}_{0}}}\left(1+\theta(\phi)\frac{V}{M_{\mathrm{pl}}^{4}}\right)$
(12)

at the time when the comoving wave number exits the horizon, where
$\mathcal{P}_{\mathcal{R}_{0}}=\frac{V^{3}}{12\pi^{2}M_{\mathrm{pl}}^{6}V_{,\phi}^{2}}$
is the power spectrum of the curvature perturbations in the minimal coupling
case. The scalar spectral index and the tensor-to-scalar ratio are given,
respectively, by S.Tsujikawa2012

$\displaystyle n_{s}\simeq
1-\frac{1}{\mathcal{A}}\left[2\epsilon_{V}\left(4-\frac{1}{\mathcal{A}}\right)-2\eta_{V}\right]\;,$
(13) $\displaystyle r\simeq\frac{16\epsilon_{V}}{\mathcal{A}}\;,$ (14)

where
$\epsilon_{V}=\frac{M_{\mathrm{pl}}^{2}}{2}\big{(}\frac{V_{,\phi}}{V}\big{)}^{2}$,
$\eta_{V}={M_{\mathrm{pl}}^{2}}\frac{V_{,\phi\phi}}{V}$ and
$\mathcal{A}=1+\frac{3}{M_{\mathrm{pl}}^{2}}\theta(\phi)H^{2}$.

For the potential of the inflaton field, we choose the simple monomial
potential

$\displaystyle V(\phi)=\lambda M_{\mathrm{pl}}^{4-p}|\phi|^{p},$ (15)

where $\lambda$ is a free parameter and the fractional power $p$ is set to be
$p=2/5$ Silverstein2008 . To amplify the curvature perturbations at the small
scales to generate a sizable amount of PBHs and at the same time to satisfy
the strong constraint on the tensor-to-scalar ratio ($r<0.036$) given by the
BICEP/Keck collaboration Ade2021 , the coupling function $\theta(\phi)$ is
assumed to take the following form Fuchengjie2019

$\displaystyle\theta(\phi)=m+\frac{\omega}{\sqrt{\kappa^{2}\left(\frac{\phi-\phi_{c}}{\sigma_{s}}\right)^{2}+1}},$
(16)

where $m$ is a coupling constant, which is introduced to reduce the tensor-to-
scalar ratio so as to be consistent with the BICEP/Keck CMB observations,
$\omega$ and $\phi_{c}$ correspond to the peak height and position of the
power spectrum of the curvature perturbations, and $\sigma_{s}$ describes the
smoothing scale around $\phi=\phi_{c}$.

At the beginning of inflation, the effect of the non-minimal derivative
coupling can be neglected since $\phi$ deviates greatly from $\phi_{c}$, and
thus the inflationary prediction corresponds to that of the standard single-
field slow-roll inflation with the simple monomial potential. The friction
will play a more and more important role with the inflaton field rolling
toward $\phi_{c}$. The large friction reduces the rolling speed of the
inflaton and leads to a period of ultra-slow-roll inflation. Since the second
term in parentheses of the r.h.s of Eq. (12) will become dominant, the power
spectrum will be enhanced. The amplitude of the power spectrum of the
curvature perturbations can be amplified to be the order of
$\mathcal{O}\left({10^{-2}}\right)$ during the ultra-slow-roll inflation. When
these enhanced curvature perturbations reenter the horizon during radiation-
or matter-dominated era, a sizable amount of PBHs will be generated.

In order to use the PBHs to explain the binary black hole events detected by
the LIGO/Virgo collaboration and the ultrashort-timescale microlensing events
in the OGLE data, and to make up all dark matter, we focus on the PBHs with
mass around $\mathcal{O}(10)M_{\odot}$, $\mathcal{O}(10^{-5})M_{\odot}$, and
$\mathcal{O}(10^{-12})M_{\odot}$, and consider three different parameter sets,
which are shown in Tab. 1. From this table, one can see that at the CMB scale
the inflationary predictions are compatible with the BICEP/Keck CMB
observations Ade2021 . And the amplitude of the power spectrum of the
curvature perturbations can be enhanced to be the
$\mathcal{O}\left({10^{-2}}\right)$ order at the small scale to generate the
abundant PBHs, as shown in Tab. 2.

##  III primordial non-gaussianity

To study the primordial non-Gaussianity of the curvature perturbations, we
need to calculate the value of the bispectrum $B_{\mathcal{R}}$, which is
related to the three-point correlation function of the curvature perturbations
C.T.Byrnes2010 ; P.Ade2016

$\displaystyle\left\langle\hat{\mathcal{R}}_{\bm{k}_{1}}\hat{\mathcal{R}}_{\bm{k}_{2}}\hat{\mathcal{R}}_{\bm{k}_{3}}\right\rangle=(2\pi)^{3}\delta^{3}\left(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3}\right)B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)\;.$
(17)

Using the in-in formula, we can calculate this three-point correlation
function and obtain the expression of the bispectrum
$B_{\mathcal{R}}(k_{1},k_{2},k_{3})$ Maldacena2003 ; FredericoArroja2011 ;
X.Chen2007 ; DavidSeery2005 ; F.Zhang2022

$\displaystyle
B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)=\Im\left[\mathcal{R}_{k_{1}}\left(t_{e}\right)\mathcal{R}_{k_{2}}\left(t_{e}\right)\mathcal{R}_{k_{3}}\left(t_{e}\right)\sum_{i=1}^{10}\mathcal{B}_{\mathcal{R}}^{i}\left(k_{1},k_{2},k_{3}\right)\right]\;.$
(18)

Here $\Im$ represents taking the imaginary part, $t_{e}$ denotes the time of
the end of inflation, and the expressions of
$\mathcal{B}_{\mathcal{R}}^{i}\left(k_{1},k_{2},k_{3}\right)$ are given in the
appendix A. Then, we can derive the non-Gaussianity parameter
$f_{\mathrm{NL}}$ C.T.Byrnes2010 ; Creminelli2002

$\displaystyle
f_{\mathrm{NL}}\left(k_{1},k_{2},k_{3}\right)=\frac{5}{6}\frac{B_{\mathcal{R}}\left(k_{1},k_{2},k_{3}\right)}{P_{\mathcal{R}}\left(k_{1}\right)P_{\mathcal{R}}\left(k_{2}\right)+P_{\mathcal{R}}\left(k_{2}\right)P_{\mathcal{R}}\left(k_{3}\right)+P_{\mathcal{R}}\left(k_{3}\right)P_{\mathcal{R}}\left(k_{1}\right)}\;,$
(19)

where $P_{\mathcal{R}}(k)=\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\mathcal{R}}(k)$.

We use the numerical method to calculate the value of
$B_{\mathcal{R}}(k_{1},k_{2},k_{3})$. Since the curvature perturbation
$\mathcal{R}$ oscillates rapidly when it is in the horizon, a cutoff
$\mathrm{e}^{\lambda k_{m}\left(\tau-\tau_{0}\right)}$ is introduced to reduce
the error in numerical calculations X.Chen2007 ; X.Chen2008 , where $k_{m}$ is
the largest value of $(k_{1},k_{2},k_{3})$, $\lambda$ determines how much the
integral will be suppressed, and $\tau_{0}$ is about several e-folding time
before the $k_{m}$ mode crosses the Hubble horizon. As the non-Gaussianity
satisfies, in the squeezed limit, the consistency relation Maldacena2003

$\displaystyle\lim_{k_{3}\rightarrow
0}f_{\mathrm{NL}}\left(k_{1},k_{2},k_{3}\right)=\frac{5}{12}(1-n_{\mathrm{s}})\quad\text{
for }k_{1}=k_{2}\gg k_{3}\;,$ (20)

it can be used to verify the accuracy of the numerical calculation.

Fig. 1 shows our numerical results in the squeezed limit for the case 1. The
solid, dashed and dotted lines represent $f_{\mathrm{NL}}$,
$\frac{5}{12}(1-n_{\mathrm{s}})$ and the power spectrum, respectively. One can
see clearly that $f_{\mathrm{NL}}$ satisfies the non-Gaussianity consistency
relation, which demonstrates fully that our numerical calculation is very
reliable. However, we must point out here that this consistency relation could
be violated if the ultra-slow-roll phase results from a flattened potential
Cai2018 ; VicenteAtal2018 ; VicenteAtal2019 . From Fig. 1, we find that at the
large scales, the power spectrum of the curvature perturbations is nearly
scale invariant, the value of $f_{\mathrm{NL}}$ is around zero, and thus non-
Gaussianity is negligible, which is the prediction of the standard slow-roll
inflation. With the decrease of scale, the power spectrum becomes to grow due
to the slowing down of the inflaton as a result of the gravitationally
enhanced friction, and accordingly the value of $f_{\mathrm{NL}}$ drops
sharply and then stabilizes at about $-0.86$. After several e-folding number,
$f_{\mathrm{NL}}$ begins to increase rapidly and reaches near zero when the
power spectrum reaches its maximum value at $k=k_{\mathrm{peak}}$. Then,
although the power spectrum decreases with the increase of $k$,
$f_{\mathrm{NL}}$ will reach its maximum value, which is about $0.47$.
Finally, with the ending of the ultra-slow-roll inflation, $f_{\mathrm{NL}}$
returns to about zero. In figure 2, we plot the $f_{\mathrm{NL}}$ in the case
of the equilateral limit ($k_{1}=k_{2}=k_{3}$), and find that it has features
similar to that of the squeezed limit case.

Figure 1:  $f_{\mathrm{NL}}$ in the squeezed limit with
$k=k_{1}=k_{2}=10^{6}k_{3}$ for the case 1. Dashed line shows
$\frac{5}{12}(1-n_{s})$ and dotted line represents the power spectrum of
curvature perturbations.  Figure 2:  $f_{\mathrm{NL}}$ in the equilateral
limit ($k=k_{1}=k_{2}=k_{3}$) for the case 1.

##  IV Effect of non-Guassianity on PBHs and SIGWS

###  IV.1 Non-Gaussian correction to the PBH abundance

When the large enough curvature perturbations reenter the Hubble horizon
during the radiation-dominated era, the gravity of overdense regions can
overcome the radiation pressure and thus these regions will collapse to form
PBHs soon after their horizon entry. The PBH mass relates with the horizon
mass at the horizon entry of perturbations with the wave number $k$:

$\displaystyle M(k)=\gamma\frac{4\pi M_{\mathrm{pl}}^{2}}{H}\simeq
M_{\odot}\left(\frac{\gamma}{0.2}\right)\left(\frac{g_{*}}{10.75}\right)^{-\frac{1}{6}}\left(\frac{k}{1.9\times
10^{6}~{}\mathrm{Mpc}^{-1}}\right)^{-2}\;,$ (21)

where $\gamma$ denotes the ratio of the PBH mass to the horizon mass and
indicates the efficiency of collapse, which is set to be
$\gamma\simeq(1/\sqrt{3})^{3}$ Carr1975  in our analysis, and $g_{*}$ is the
effective degrees of freedom in the energy densities at the PBH formation. We
adopt $g_{*}=106.75$ since the PBHs are assumed to form deep in the radiation-
dominated.

Based on the Press-Schechter theory Tada2019 ; Young2014 , the production rate
of PBHs with mass $M(k)$ has the form

$\displaystyle\beta(M)=\int_{\delta_{c}}\frac{d\delta}{\sqrt{2\pi\sigma^{2}(M)}}e^{-\frac{\delta^{2}}{2\sigma^{2}(M)}}=\frac{1}{2}\operatorname{erfc}\left(\frac{\delta_{c}}{\sqrt{2\sigma^{2}(M)}}\right),$
(22)

after assuming that the probability distribution function of perturbations is
Gaussian, where $\mathrm{erfc}$ is the complementary error function.
$\delta_{c}$ is the threshold of the density perturbations for the PBH
formation, which is chosen to be $\delta_{c}\simeq 0.4$ Musco2013 ; Harada2013
in our calculation of PBHs abundance. $\sigma^{2}(M)$ has the form

$\displaystyle\sigma^{2}(M(k))=\frac{16}{81}\int d\ln
q\;W^{2}\left(qk^{-1}\right)\left(qk^{-1}\right)^{4}\mathcal{P}_{\mathcal{R}}(q)\;,$
(23)

which represents the coarse-grained density contrast with the smoothing scale
$k$. Here $W$ is the window function. The current fractional energy density of
PBHs with mass $M$ in dark matter is

$\displaystyle
f(M)\equiv\frac{1}{\Omega_{\mathrm{DM}}}\frac{d\Omega_{\mathrm{PBH}}}{d\ln
M}\simeq\frac{\beta(M)}{1.84\times
10^{-8}}\left(\frac{\gamma}{0.2}\right)^{\frac{3}{2}}\left(\frac{10.75}{g_{*}}\right)^{\frac{1}{4}}\left(\frac{0.12}{\Omega_{\mathrm{DM}}h^{2}}\right)\left(\frac{M}{M_{\odot}}\right)^{-\frac{1}{2}},$
(24)

where $\Omega_{\mathrm{DM}}$ is the current density parameter of dark matter,
which is given to be $\Omega_{\mathrm{DM}}h^{2}\simeq 0.12$ by the Planck 2018
observations Aghanim2020 . For the Gaussian distribution of the curvature
perturbations, we obtain PBHs with masses around
$\mathcal{O}\left({10}\right)M_{\odot}$,
$\mathcal{O}\left({10^{-6}}\right)M_{\odot}$, and
$\mathcal{O}\left({10^{-13}}\right)M_{\odot}$, respectively and their
corresponding abundances, which are shown in Tab. 2.

When the effect of non-Guassianity of the curvature perturbations on the PBH
abundance is considered, the mass fraction $\beta$ is corrected to be
G.Franciolini2018 ; Riccardi2021

$\displaystyle\beta=\mathrm{e}^{\Delta_{3}}\beta^{G}\;,$ (25)

where $\Delta_{3}$ is the 3rd cumulant, which has the form

$\displaystyle\Delta_{3}=\frac{1}{3!}\left(\frac{\delta_{c}}{\sigma}\right)^{2}S_{3}\delta_{c},$
(26)

with $S_{3}$ being

$\displaystyle
S_{3}=\frac{\left\langle\delta_{R}(\bm{x})\delta_{R}(\bm{x})\delta_{R}(\bm{x})\right\rangle}{\sigma^{4}}\;.$
(27)

For the Gaussian window function,
$\left\langle\delta_{R}(\bm{x})\delta_{R}(\bm{x})\delta_{R}(\bm{x})\right\rangle$
can be obtained through calculating

$\displaystyle\left\langle\delta_{R}(\bm{x})\delta_{R}(\bm{x})\delta_{R}(\bm{x})\right\rangle$
$\displaystyle=$
$\displaystyle-64\left(\frac{4}{9}\right)^{3}\frac{2}{(2\pi)^{4}}k^{6}$
$\displaystyle\times$
$\displaystyle\int_{0}^{\infty}du\int_{0}^{\infty}dv\int_{|u-v|}^{u+v}dwu^{3}v^{3}w^{3}\mathrm{e}^{-u^{2}}\mathrm{e}^{-v^{2}}\mathrm{e}^{-w^{2}}B_{\mathcal{R}}(\sqrt{2}uk,\sqrt{2}vk,\sqrt{2}wk).$

In Zhang2021 , it has been found that $\Delta_{3}$ can be approximated as

$\displaystyle\Delta_{3}\approx\Delta_{3}^{a}(k_{\text{peak}})=23\frac{\delta_{c}^{3}}{\mathcal{P}_{\mathcal{R}}\left(k_{\mathrm{peak}}\right)}~{}f_{\text{NL
}}(k_{\mathrm{peak}},k_{\mathrm{peak}},k_{\mathrm{peak}}).$ (29)

We numerically compute $\Delta_{3}$ and its approximation $\Delta_{3}^{a}$.
The results are shown in Tab. 2 for three different cases. It is easy to see
that $\Delta_{3}$ is of the order $\mathcal{O}(1)$, which is much larger than
$\Delta_{3}^{a}$. Thus, the approximation given in Zhang2021  is invalid for
the model considered in the present paper. We find that for all the cases the
non-Gaussianity promotes the formation of PBHs since their 3rd cumulants
$\Delta_{3}$ are positive. The effect of non-Gaussianity on the PBH abundance
is non-negligible since the value of $\beta$ is much larger than $\beta^{G}$.

Table 2: The numerical results for three cases given in Tab. 1.  | $k_{\mathrm{peak}}/\mathrm{Mpc}^{-1}$ |  $\mathcal{P}_{\mathcal{R}}/10^{-2}$ |  $f_{\mathrm{PBH}}^{G}$ |  $M_{\mathrm{PBH}}/M_{\odot}$ | $f_{\mathrm{NL}}$ | $\Delta_{3}^{a}$ | $\Delta_{3}$
---|---|---|---|---|---|---|---
Case 1 | $4.46\times 10^{5}$ | $4.91$ | $0.00455$ | $16.57$ | $~{}~{}0.0205$ | $~{}~{}0.619$ | $~{}4.218$
Case 2 | $4.31\times 10^{8}$ | $3.87$ | $0.0178$ | $~{}~{}1.91\times 10^{-5}$ | $~{}~{}0.0234$ | $~{}~{}0.890$ | $~{}4.965$
Case 3 | $1.77\times 10^{12}$ | $3.18$ | $0.924$ |  $1.12\times 10^{-12}$ | $-0.0052$ | $-0.240$ | $~{}6.297$

###  IV.2 Non-Gaussian correction to energy density of SIGWs

Associated with the PHB formation, the enhanced curvature perturbations will
lead to the large scalar metric perturbations, which become the significant GW
source and emit abundant SIGWs. The current energy spectra of SIGWs can be
expressed as Kohri2018 ; Inomata2019

$\displaystyle\Omega_{\mathrm{GW},0}h^{2}=0.83\left(\frac{g_{*}}{10.75}\right)^{-1/3}\Omega_{\mathrm{r},0}h^{2}\Omega_{\mathrm{GW}}\left(\tau_{c},k\right),$
(30)

where $\Omega_{\mathrm{r},0}h^{2}\simeq 4.2\times 10^{5}$ is the current
density parameter of radiation, and

$\displaystyle\Omega_{\mathrm{GW}}\left(\tau_{c},k\right)=$
$\displaystyle\frac{1}{12}\int_{0}^{\infty}dv\int_{|1-v|}^{|1+v|}du\left(\frac{4v^{2}-\left(1+v^{2}-u^{2}\right)^{2}}{4uv}\right)^{2}\mathcal{P}_{\mathcal{R}}(ku)\mathcal{P}_{\mathcal{R}}(kv)$
$\displaystyle\times\left(\frac{3}{4u^{3}v^{3}}\right)^{2}\left(u^{2}+v^{2}-3\right)^{2}$
$\displaystyle\times\bigg{\\{}\left[-4uv+\left(u^{2}+v^{2}-3\right)\ln\left|\frac{3-(u+v)^{2}}{3-(u-v)^{2}}\right|\right]^{2}$
$\displaystyle+\pi^{2}\left(u^{2}+v^{2}-3\right)^{2}\Theta(v+u-\sqrt{3})\bigg{\\}}\,,$
(31)

where $\tau_{c}$ represents the time when $\Omega_{\mathrm{GW}}$ stops to grow
and $\Theta$ is the Heaviside theta function.

When the non-Gaussianity of the curvature perturbations is considered,
$\mathcal{R}(\bm{x})$ has the expression Verde2000 ; Komatsu2001

$\displaystyle\mathcal{R}(\bm{x})=\mathcal{R}^{G}(\bm{x})+\frac{3}{5}f_{\mathrm{NL}}\left(\mathcal{R}^{G}(\bm{x})^{2}-\left\langle\mathcal{R}^{G}(\bm{x})^{2}\right\rangle\right)\,.$
(32)

Clearly the curvature perturbation $\mathcal{R}$ consists of a Gaussian part
$\mathcal{R}^{G}$ and a non-Gaussian one. When this non-Gaussian correction is
included, the power spectrum of the curvature perturbations should be modified
to be

$\displaystyle\mathcal{P}_{\mathcal{R}}(k)=\mathcal{P}_{\mathcal{R}}^{G}(k)+\mathcal{P}_{\mathcal{R}}^{NG}(k),$
(33)

where

$\displaystyle\mathcal{P}_{\mathcal{R}}^{NG}(k)=\left(\frac{3}{5}\right)^{2}\frac{k^{3}}{2\pi}f_{\mathrm{NL}}^{2}\int
d^{3}\bm{p}\frac{\mathcal{P}_{\mathcal{R}}^{G}(p)}{p^{3}}\frac{\mathcal{P}_{\mathcal{R}}^{G}(|\bm{k}-\bm{p}|)}{|\bm{k}-\bm{p}|^{3}}.$
(34)

For the model considered in this paper, the result in the preceding subsection
has shown that the absolute value of $f_{\mathrm{NL}}$ is less than one and it
is near zero when $k=k_{\mathrm{peak}}$. Furthermore the maximum value of
$\mathcal{P}_{\mathcal{R}}^{G}$ has the order of $\mathcal{O}(0.01)$. Thus, we
can assess easily that $\mathcal{P}_{\mathcal{R}}^{NG}(k)$ is much less than
$\mathcal{P}_{\mathcal{R}}^{G}(k)$ since the order of its maximum should be
less than $\mathcal{O}(10^{-4})$, which indicates that the contribution of
non-Gaussianity of the curvature perturbations on the energy density of SIGWs
is negligible.

##  V conclusions

To generate a sizable amount of PBHs requires that the amplitude of the power
spectrum of the curvature perturbations is enhanced to reach the
$\mathcal{O}(0.01)$ order. A simple way to enhance the curvature perturbations
during inflation is to reduce the rolling speed of inflaton to achieve an
ultra-slow-roll inflation, which can be realized by flattening the
inflationary potential or increasing the gravitational friction. Since the
ultra-slow-roll inflation deviates apparently from the standard slow-roll
inflation, the non-Gaussianity of the curvature perturbations might be very
large and has significant effects on the abundance of PBHs and the energy
density of SIGWs although it is negligible in the standard slow-roll
inflation.

In this paper we study the non-Gaussianity of the curvature perturbations in
the ultra-slow-roll inflation resulting from gravitationally enhanced
friction. We find that at the large scales where the power spectrum of the
curvature perturbations is nearly scale invariant, the non-Gaussianity is
negligible. The power spectrum grows with the decrease of scale due to that
the friction slows down the inflaton, and correspondingly the value of the
non-Gaussianity parameter $f_{\mathrm{NL}}$ drops sharply and then stabilizes
at a value in several e-folding number. Before the power spectrum reaches its
peak, $f_{\mathrm{NL}}$ begins to increase rapidly. When the power spectrum is
at its maximum value, we find that $f_{\mathrm{NL}}$ is nearly zero. Then,
$f_{\mathrm{NL}}$ will reach its maximum value with the increase of wave
number $k$. Finally, with the ending of the ultra-slow-roll inflation,
$f_{\mathrm{NL}}$ returns to about zero. For three different cases, which
correspond to that the PBHs can be used to explain the LIGO/Virgo GW events
and the six ultrashort-timescale microlensing events in the OGLE data, and
make up all dark matter, respectively, we obtain that the non-Gaussianity will
promote the generation of PBHs, while its influence on SIWGs is negligible.

##  Appendix A The expressions of
$\mathcal{B}_{\mathcal{R}}^{i}\left(k_{1},k_{2},k_{3}\right)$ in Eq. (18)

In order to compute bispectrum, we need derive the cubic action of the
curvature perturbations from the action given in Eq. (1) A.D.Felice2011 ;
Felice2011

$\displaystyle S_{3}$ $\displaystyle=$ $\displaystyle\int
dtd^{3}x\left\\{a^{3}\mathcal{C}_{1}\mathcal{R}\dot{\mathcal{R}}^{2}+a\mathcal{C}_{2}\mathcal{R}(\partial\mathcal{R})^{2}+a^{3}\mathcal{C}_{3}\dot{\mathcal{R}}^{3}+a^{3}\mathcal{C}_{4}\dot{\mathcal{R}}\left(\partial_{i}\mathcal{R}\right)\left(\partial_{i}\chi\right)\right.$
(35)
$\displaystyle+a^{3}\mathcal{C}_{5}\partial^{2}\mathcal{R}(\partial\chi)^{2}+a\mathcal{C}_{6}\dot{\mathcal{R}}^{2}\partial^{2}\mathcal{R}+\left(\mathcal{C}_{7}/a\right)\left[\partial^{2}\mathcal{R}(\partial\mathcal{R})^{2}-\mathcal{R}\partial_{i}\partial_{j}\left(\partial_{i}\mathcal{R}\right)\left(\partial_{j}\mathcal{R}\right)\right]$
$\displaystyle+a\mathcal{C}_{8}\left[\partial^{2}\mathcal{R}\partial_{i}\mathcal{R}\partial_{i}\chi-\mathcal{R}\partial_{i}\partial_{j}\left(\partial_{i}\mathcal{R}\right)\left(\partial_{j}\chi\right)\right]+\left.\mathcal{F}_{1}\frac{\delta\mathcal{L}_{2}}{\delta\mathcal{R}}\right\\},$

where

$\displaystyle\mathcal{F}_{1}$ $\displaystyle=$
$\displaystyle\frac{A_{4}}{w_{1}^{2}}\left\\{\left(\partial_{k}\mathcal{R}\right)\left(\partial_{k}\chi\right)-\partial^{-2}\partial_{i}\partial_{j}\left[\left(\partial_{i}\mathcal{R}\right)\left(\partial_{j}\chi\right)\right]\right\\}$
$\displaystyle+q_{1}\mathcal{R}\dot{\mathcal{R}}-\frac{A_{4}}{a^{2}w_{2}}\times\left\\{(\partial\mathcal{R})^{2}-\partial^{-2}\partial_{i}\partial_{j}\left[\left(\partial_{i}\mathcal{R}\right)\left(\partial_{j}\mathcal{R}\right)\right]\right\\},$

$\partial^{2}\chi=Q\dot{\mathcal{R}}$, $\mathcal{L}_{2}$ is quadratic
Lagrangian given in Eq. (7), $w_{i}$ and $Q$ are given in Eqs. (II) and (8),
respectively, and the dimensionless coefficients $\mathcal{C}_{i}$ with
$i=1-8$ are

$\displaystyle\mathcal{C}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{M_{\mathrm{pl}}^{2}}\left[3Q+q_{1}(\dot{Q}+3HQ)-Q\dot{q}_{1}\right],$
(36) $\displaystyle\mathcal{C}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{M_{\mathrm{pl}}^{2}}\left[A_{5}+\frac{1}{a}\frac{d}{dt}\left(\frac{2aQw_{1}}{w_{2}}\right)\right],$
(37) $\displaystyle\mathcal{C}_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{M_{\mathrm{pl}}}\left[A_{1}+A_{3}\frac{Q}{w_{1}}-q_{1}Q\right],$
(38) $\displaystyle\mathcal{C}_{4}$ $\displaystyle=$
$\displaystyle\frac{Q}{w_{1}}\left[-\frac{1}{2}-w_{1}\frac{d}{dt}\left(\frac{A_{4}}{w_{1}^{2}}\right)+\frac{3HA_{4}}{w_{1}}\right],$
(39) $\displaystyle\mathcal{C}_{5}$ $\displaystyle=$
$\displaystyle\frac{M_{\mathrm{pl}}^{2}}{2}\left[\frac{3}{2w_{1}}-\frac{d}{dt}\left(\frac{A_{4}}{w_{1}^{2}}\right)+\frac{3HA_{4}}{w_{1}^{2}}\right],$
(40) $\displaystyle\mathcal{C}_{6}$ $\displaystyle=$ $\displaystyle
A_{2}-\frac{2w_{1}A_{3}}{w_{2}},$ (41) $\displaystyle\mathcal{C}_{7}$
$\displaystyle=$ $\displaystyle q_{3}+\frac{2A_{4}Qc_{s}^{2}}{w_{2}},$ (42)
$\displaystyle\mathcal{C}_{8}$ $\displaystyle=$ $\displaystyle
M_{\mathrm{pl}}\left[\frac{q_{2}}{2}-\frac{2c_{s}^{2}A_{4}Q}{w_{1}^{2}}\right].$
(43)

Here

$\displaystyle A_{1}$ $\displaystyle=$
$\displaystyle\frac{3w_{1}}{w^{3}_{2}}\bigg{[}8H^{2}\left(M_{\mathrm{pl}}^{2}-5\dot{\phi}^{2}\theta\left(\phi\right)\right)w_{1}^{2}+8H\left(3\dot{\phi}^{2}\theta\left(\phi\right)-M_{\mathrm{pl}}^{2}\right)w_{1}w_{2}+\left(2M_{\mathrm{pl}}^{2}-3\dot{\phi}^{2}\theta\left(\phi\right)\right)w_{2}^{2}\bigg{]},$
$\displaystyle A_{2}$ $\displaystyle=$
$\displaystyle-4\dot{\phi}^{2}\theta(\phi)\frac{w_{1}^{2}}{w^{2}_{2}},$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\frac{2w_{1}}{w^{2}_{2}}\left[4Hw_{1}\left(M_{\mathrm{pl}}^{2}-3\dot{\phi}^{2}\theta\left(\phi\right)\right)+w_{2}\left(3\dot{\phi}^{2}\theta\left(\phi\right)-2M_{\mathrm{pl}}^{2}\right)\right],$
$\displaystyle A_{4}$ $\displaystyle=$
$\displaystyle\frac{w_{1}}{2w_{2}}\left(3\dot{\phi}^{2}\theta\left(\phi\right)-2M_{\mathrm{pl}}^{2}\right),$
$\displaystyle A_{5}$ $\displaystyle=$
$\displaystyle\frac{2\dot{w}_{2}w_{1}^{2}+w_{2}\left(w_{2}w_{4}-4\dot{w}_{1}w_{1}-2Hw_{1}^{2}\right)}{w_{2}^{2}},$

and

$\displaystyle q_{1}$ $\displaystyle=$
$\displaystyle-\frac{2w_{1}}{c_{s}^{2}w_{2}},$ $\displaystyle q_{2}$
$\displaystyle=$ $\displaystyle
a^{2}\frac{d}{dt}\left[\frac{w_{1}\left(4M_{\mathrm{pl}}^{2}-6\dot{\phi}^{2}\theta\left(\phi\right)\right)}{a^{2}w_{2}^{2}}\right]-\frac{4w_{1}}{w_{2}},$
$\displaystyle q_{3}$ $\displaystyle=$
$\displaystyle\frac{2w_{1}^{3}}{3w_{2}^{2}}-\frac{a}{3}\frac{d}{dt}\left[\frac{2w_{1}^{3}\left(3\dot{\phi}^{2}\theta(\phi)-2M_{\mathrm{pl}}^{2}\right)}{aw_{2}^{3}}\right].$

Using the in-in formula, one can obtain the three-point correlation function
from the cubic action of the curvature perturbations. The analytical
expression is shown in Eq. (18), in which
$\mathcal{B}_{\mathcal{R}}^{i}\left(k_{1},k_{2},k_{3}\right)$ have the forms

$\displaystyle\mathcal{B}_{\mathcal{R}}^{1}\left(k_{1},k_{2},k_{3}\right)=-4\int_{t_{i}}^{t_{e}}dt\;a^{3}\;\mathcal{C}_{1}\left(\mathcal{R}_{k_{1}}^{*}(t)\dot{\mathcal{R}}_{k_{2}}^{*}(t)\dot{\mathcal{R}}^{*}_{k_{3}}(t)+\text{
perm }\right),$ (44)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{2}\left(k_{1},k_{2},k_{3}\right)=4\int_{t_{i}}^{t_{e}}dt\;a\;\mathcal{C}_{2}\left[\left(\bm{k}_{1}\cdot\bm{k}_{2}+\bm{k}_{1}\cdot\bm{k}_{3}+\bm{k}_{2}\cdot\bm{k}_{3}\right)\mathcal{R}_{k_{1}}^{*}(t)\mathcal{R}_{k_{2}}^{*}(t)\mathcal{R}_{k_{3}}^{*}(t)\right],$
(45)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{3}\left(k_{1},k_{2},k_{3}\right)=-12\int_{t_{i}}^{t_{e}}dt\;a^{3}\;\mathcal{C}_{3}\dot{\mathcal{R}}_{k_{1}}^{*}(t)\dot{\mathcal{R}}_{k_{2}}^{*}(t)\dot{\mathcal{R}}_{k_{3}}^{*}(t),$
(46)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{4}\left(k_{1},k_{2},k_{3}\right)=-2\int_{t_{i}}^{t_{e}}dt\;a^{3}\;\mathcal{C}_{4}\;Q\left[\left(\frac{\bm{k}_{1}\cdot\bm{k}_{2}}{k_{2}^{2}}+\frac{\bm{k}_{1}\cdot\bm{k}_{3}}{k_{3}^{2}}\right)\mathcal{R}_{k_{1}}^{*}(t)\dot{\mathcal{R}}_{k_{2}}^{*}(t)\dot{\mathcal{R}}_{k_{3}}^{*}(t)+\text{
perm }\right],$ (47)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{5}\left(k_{1},k_{2},k_{3}\right)=-4\int_{t_{i}}^{t_{e}}dt\;a^{3}\;\mathcal{C}_{5}\;Q^{2}\left[\frac{k_{1}^{2}\bm{k}_{2}\cdot\bm{k}_{3}}{k_{2}^{2}k_{3}^{2}}\mathcal{R}_{k_{1}}^{*}(t)\dot{\mathcal{R}}_{k_{2}}^{*}(t)\dot{\mathcal{R}}_{k_{3}}^{*}(t)+\text{
perm }\right],$ (48)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{6}\left(k_{1},k_{2},k_{3}\right)=4\int_{t_{i}}^{t_{e}}dt\;a\;\mathcal{C}_{6}\left[k_{1}^{2}\mathcal{R}_{k_{1}}^{*}(t)\dot{\mathcal{R}}^{*}k_{2}(t)\dot{\mathcal{R}}^{*}k_{3}(t)+\text{
perm }\right],$ (49)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{7}\left(k_{1},k_{2},k_{3}\right)=-4\int_{t_{i}}^{t_{e}}dt\;\mathcal{C}_{7}/a\left[\left(k_{1}^{2}\bm{k}_{2}\cdot\bm{k}_{3}+k_{2}^{2}\bm{k}_{1}\cdot\bm{k}_{3}+k_{3}^{2}\bm{k}_{1}\cdot\bm{k}_{2}\right)\mathcal{R}_{k_{1}}^{*}(t)\mathcal{R}_{k_{2}}^{*}(t)\mathcal{R}_{k_{3}}^{*}(t)\right],$
(50)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{8}\left(k_{1},k_{2},k_{3}\right)=2\int_{t_{i}}^{t_{e}}dt\;\mathcal{C}_{7}/a$
$\displaystyle{\left[\left(k_{2}^{2}\bm{k}_{2}\cdot\bm{k}_{3}+k_{3}^{2}\bm{k}_{3}\cdot\bm{k}_{2}+k_{1}^{2}\bm{k}_{1}\cdot\bm{k}_{2}+k_{1}^{2}\bm{k}_{1}\cdot\bm{k}_{3}\right.\right.}$
$\displaystyle\left.\left.+k_{2}^{2}\bm{k}_{2}\cdot\bm{k}_{1}+k_{3}^{2}\bm{k}_{3}\cdot\bm{k}_{1}\right)\mathcal{R}_{k_{1}}^{*}(t)\mathcal{R}_{k_{2}}^{*}(t)\mathcal{R}_{k_{3}}^{*}(t)\right],$
(51)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{9}\left(k_{1},k_{2},k_{3}\right)=2\int_{t_{i}}^{t_{e}}dt\;a\;\mathcal{C}_{8}\;Q\left[\left(\frac{k_{1}^{2}\bm{k}_{2}\cdot\bm{k}_{3}}{k_{3}^{2}}+\frac{k_{2}^{2}\bm{k}_{1}\cdot\bm{k}_{3}}{k_{3}^{2}}\right)\mathcal{R}_{k_{1}}^{*}(t)\mathcal{R}_{k_{2}}^{*}(t)\dot{\mathcal{R}}_{k_{3}}^{*}(t)+\text{
perm }\right],$ (52)
$\displaystyle\mathcal{B}_{\mathcal{R}}^{10}\left(k_{1},k_{2},k_{3}\right)=-2\int_{t_{i}}^{t_{e}}dt\;a\;\mathcal{C}_{8}\;Q\left[\left(\frac{k_{2}^{2}\bm{k}_{2}\cdot\bm{k}_{3}}{k_{3}^{2}}+\frac{k_{1}^{2}\bm{k}_{1}\cdot\bm{k}_{3}}{k_{3}^{2}}\right)\mathcal{R}_{k_{1}}^{*}(t)\mathcal{R}_{k_{2}}^{*}(t)\dot{\mathcal{R}}_{k_{3}}^{*}(t)+\operatorname{perm}\right].$
(53)

###### Acknowledgements.

We appreciate very much the insightful comments and helpful suggestions by the
anonymous referee. L. Chen is grateful to Dr. Fengge Zhang for his valuable
help in the non-Gaussian numerical computation. This work is supported by the
National Key Research and Development Program of China Grant No.
2020YFC2201502, and by the National Natural Science Foundation of China under
Grants No. 12275080, No. 12075084, and No. 11805063.

## References

  * (1) A. H. Guth, Phys. Rev. D 23, 347 (1981).
  * (2) A. D. Linde, Phys. Lett. 108B, 389 (1982).
  * (3) A. A. Starobinsky, Phys. Lett. 91B, 99 (1980).
  * (4) A. Albrecht and P. J. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982).
  * (5) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020).
  * (6) S. Hawking, Mon. Not. R. Astron. Soc. 152, 75 (1971).
  * (7) B. J. Carr and S. W. Hawking, Mon. Not. R. Astron. Soc. 168, 399 (1974).
  * (8) B. J. Carr,  Astrophys. J. 201, 1 (1975).
  * (9) M. Y. Khlopov, Res. Astron. Astrophys. 10, 495 (2010).
  * (10) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016).
  * (11) S. Clesse and J. García-Bellido, Phys. Dark Univ. 15, 142 (2017).
  * (12) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016).
  * (13) B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016).
  * (14) P. Mróz, A. Udalski, J. Skowron, et al., Nature (London) 548, 183 (2017).
  * (15) H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Phys. Rev. D 99, 083503 (2019).
  * (16) A. Katz, J. Kopp, S. Sibiryakov, and W. Xue, J. Cosmol. Astropart. Phys. 12 (2018) 005.
  * (17) A. Barnacka, J. F. Glicenstein, and R. Moderski, Phys. Rev. D 86, 043001 (2012).
  * (18) P. W. Graham, S. Rajendran, and J. Varela, Phys. Rev. D 92, 063007 (2015).
  * (19) H. Niikura, M. Takada, N. Yasuda, et al., Nat. Astron. 3, 524 (2019).
  * (20) K. M. Belotsky, A. D. Dmitriev, E. A. Esipova, V. A. Gani, A. V. Grobov, M. Y. Khlopov, A.A.Kirillov, S. G. Rubin, and I. V. Svadkovsky,  Mod. Phys. Lett. A  29, 1440005 (2014).
  * (21) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6 (2017).
  * (22) H. Motohashi and W. Hu, Phys. Rev. D 96, 063503 (2017).
  * (23) J. M. Ezquiaga, J. Garcia-Bellido, and E. Ruiz Morales, Phys. Lett. 039B, 11 (2017).
  * (24) H. Di and Y. Gong, J. Cosmol. Astropart. Phys. 07 (2018) 007.
  * (25) G. Ballesteros and M. Taoso, Phys. Rev. D 97, 023501 (2018).
  * (26) I. Dalianis, A. Kehagias, and G. Tringas, J. Cosmol. Astropart. Phys. 01 (2019) 037.
  * (27) T. J. Gao and Z. K. Guo, Phys. Rev. D 98, 063526 (2018).
  * (28) M. Drees and Y. Xu, Eur. Phys. J. C 81, 182 (2021).
  * (29) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 102, 043527 (2020).
  * (30) W. Xu, J. Liu, T. Gao, and Z. Guo, Phys. Rev. D 101, 023505 (2020).
  * (31) J. Lin, Q. Gao, Y. Gong, Y. Lu, and C. Zhang, Phys. Rev. D 101, 103515 (2020).
  * (32) I. Dalianis and K. Kritos, Phys. Rev. D 103, 023505 (2021).
  * (33) Z. Yi, Y. Gong, B. Wang, and Z. Zhu, Phys. Rev. D 103, 063535 (2021).
  * (34) Q. Gao, Y. Gong, and Z. Yi, Nucl. Phys. B 969, 115480 (2021).
  * (35) Z. Yi, Q. Gao, Y. Gong, and Z. Zhu, Phys. Rev. D 103, 063534 (2021).
  * (36) T. Gao and X. Yang, Eur. Phys. J. C 81, 494 (2021).
  * (37) M. Solbi and K. Karami, J. Cosmol. Astropart. Phys. 08 (2021) 056.
  * (38) Q. Gao, Sci. China Phys. Mech. Astron. 64, 280411 (2021).
  * (39) M. Solbi and K. Karami, Eur. Phys. J. C 81, 884 (2021).
  * (40) R. Zheng, J. Shi, and T. Qiu, Chin. Phys. C 46 045103.
  * (41) Z. Teimoori, K. Rezazadeh, M. A. Rasheed, and K. Karami, J. Cosmol. Astropart. Phys. 10 (2021) 018.
  * (42) R. Cai, C. Chen, and C. Fu, Phys. Rev. D 104, 083537 (2021).
  * (43) Q. Wang, Y. Liu, B. Su, and N. Li, Phys. Rev. D 104, 083546 (2021).
  * (44) A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen, and H. Veermäe, arXiv.2205.13540.
  * (45) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 100, 063532 (2019).
  * (46) C. Fu, P. Wu, and H. Yu, Phys. Rev. D 101, 023529 (2020).
  * (47) I. Dalianis, S. Karydas, and E. Papantonopoulos, J. Cosmol. Astropart. Phys. 06 (2020) 040.
  * (48) Z. Teimoori, K. Rezazadeh, and K. Karami, Astrophys. J. 915, 118 (2021).
  * (49) S. Heydari, and K. Karami,  Eur. Phys. J. C 82, 83 (2022).
  * (50) S. Heydari, and K. Karami, J. Cosmol. Astropart. Phys. 03 (2022) 033.
  * (51) Y. F. Cai, X. Tong, D. G. Wang and S. F. Yan, Phys. Rev. Lett. 121 (2018) 081306.
  * (52) Y. F. Cai, C. Chen, X. Tong, D. G. Wang and S. F. Yan, Phys. Rev. D 100 (2019) 043518.
  * (53) C. Chen and Y. F. Cai, J. Cosmol. Astropart. Phys. 10 (2019) 068.
  * (54) C. Chen, X. H. Ma and Y. F. Cai, Phys. Rev. D 102 (2020) 063526.
  * (55) A. Addazi, S. Capozziello, and Q. Gan, arXiv: 2204.07668
  * (56) R. Cai, Z. Guo, J. Liu, L. Liu, and X. Yang, J. Cosmol. Astropart. Phys. 06 (2020) 013
  * (57) P. Amaro-Seoane, H. Audley, S. Babak, et al. (LISA), arXiv:1702.00786.
  * (58) W. R. Hu and Y. L. Wu, Natl. Sci. Rev. 4, 685 (2017).
  * (59) J. Luo, L. S. Chen, H. Z. Duan, et al. (TianQin), Class. Quant. Grav. 33, 035010 (2016).
  * (60) R. D. Ferdman, R. van Haasteren, C. G. Bassa, et al., Class. Quant. Grav. 27, 084014 (2010).
  * (61) G. Hobbs, A. Archibald2, Z. Arzoumanian, et al., Class. Quant. Grav. 27, 084013 (2010).
  * (62) M. A. McLaughlin, Class. Quant. Grav. 30, 224008 (2013).
  * (63) G. Hobbs, Class. Quant. Grav. 30, 224007 (2013).
  * (64) R. G. Cai, S. Pi, and M. Sasaki, Phys. Rev. Lett. 122, 201101 (2019).
  * (65) F. Zhang, Phys. Rev. D 105, 063539 (2022).
  * (66) F. Zhang, J. Cosmol. Astropart. Phys. 04 (2021) 045.
  * (67) F. Zhang, J. Lin, and Y. Lu,  Phys. Rev. D 104, 063515 (2021).
  * (68) S. Passaglia, W. Hu, and H. Motohashi, Phys. Rev. D 99, 043536 (2019).
  * (69) C. M. Yoo, J. O. Gong, and S. Yokoyama, J. Cosmol. Astropart. Phys. 09 (2019) 033.
  * (70) G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, J. Cosmol. Astropart. Phys. 03 (2018) 016.
  * (71) M. W. Davies, P. Carrilho, D. J. Mulryne, J. Cosmol. Astropart. Phys. 06 (2022) 019.
  * (72) Q. G. Huang and Y. Wang, J. Cosmol. Astropart. Phys. 06 (2013) 035.
  * (73) B. Rafael, M. Sander, P. G. A, Bastian, J. Cosmol. Astropart. Phys. 05 (2018) 025.
  * (74) Y. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D. Wang, and Z. Wang, J. Cosmol. Astropart. Phys. 05 (2018) 012.
  * (75) V. Atal and C. Germani, Phys. Dark Univ. 24, 1002757 (2019).
  * (76) V. Atal, J. Garriga and A. Marcos-Caballero J. Cosmol. Astropart. Phys. 09 (2019) 073.
  * (77) C. Germani and A. Kehagias, Phys. Rev. Lett. 106, 161302 (2011).
  * (78) C. Germani and Y. Watanabe, J. Cosmol. Astropart. Phys. 07 (2011) 031.
  * (79) A. D. Felice and S.Tsujikawa, J. Cosmol. Astropart. Phys. 04 (2011) 029.
  * (80) S. Tsujikawa, Phys. Rev. D 85, 083518 (2012).
  * (81) T. Kobayashi, M. Yamaguchi, and J. Yokoyama, Prog. Theor. Phys. 126, 511 (2011).
  * (82) E. Silverstein and A. Westphal, Phys. Rev. D 78, 106003 (2008).
  * (83) P. A. R. Ade, Z. Ahmed, M. Amiri, et al., Phys. Rev. Lett. 127, 151301 (2021).
  * (84) C. T. Byrnes, M. Gerstenlauer, S. Nurmi, G. Tasinato, and D. Wands, J. Cosmol. Astropart. Phys. 10 (2010) 004.
  * (85) P. A. R. Ade, N. Aghanim, M. Arnaud, et al, Astron. Astrophys. 594, A17 (2016).
  * (86) D. Seery and J. E. Lidsey, J. Cosmol. Astropart. Phys. 06 (2005) 003.
  * (87) X. Chen, R. Easther, and E. A. Lim, Cosmol. Astropart. Phys. 06 (2007) 023.
  * (88) F. Arroja and T. Tanaka, J. Cosmol. Astropart. Phys. 05 (2011) 005.
  * (89) J. M. Maldacena, J. High Energ. Phys. 05 (2003) 013.
  * (90) P. Creminelli, L. Senatore, M. Zaldarriaga, and M. Tegmark, J. Cosmol. Astropart. Phys. 03 (2007) 005.
  * (91) X. Chen, R. Easther, and E. A. Lim, J. Cosmol. Astropart. Phys. 04 (2008) 010.
  * (92) Y. Tada and S. Yokoyama, Phys. Rev. D 100, 023537 (2019).
  * (93) S. Young, C. T. Byrnes, and M. Sasaki, J. Cosmol. Astropart. Phys. 07 (2014) 045.
  * (94) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013).
  * (95) T. Harada, C. M. Yoo, and K. Kohri, Phys. Rev. D 88, 084051 (2013).
  * (96) F. Riccardi, M. Taoso, and A. Urbano, J. Cosmol. Astropart. Phys. 08 (2021) 060.
  * (97) K. Inomata and T. Nakama, Phys. Rev. D 99, 043511 (2019).
  * (98) K. Kohri and T. Terada, Phys. Rev. D 97, 123532 (2018).
  * (99) L. Verde, L. M. Wang, A. Heavens, and M. Kamionkowski, Mon. Not. Roy. Astron. Soc. 313, L141 (2000).
  * (100) E. Komatsu and D. N. Spergel, Phys. Rev. D 63, 063002 (2001).
  * (101) A. De Felice and S. Tsujikawa, Phys. Rev. D 84, 083504 (2011).


# On some invariants of Birkhoff billiards under conjugacy

V. Kaloshin111Vadim Kaloshin: University of Maryland, College Park, MD, USA &
Institute of Science and Technology Austria, Am Campus 1, 3400 Klosterneuburg,
Austria<EMAIL_ADDRESS>, C. E. Koudjinan222Comlan Edmond
Koudjinan: Institute of Science and Technology Austria (IST Austria), Am
Campus 1, 3400 Klosterneuburg, Austria<EMAIL_ADDRESS>
###### Abstract

In the class of strictly convex smooth boundaries, each of which not having
strip around its boundary foliated by invariant curves, we prove that the
Taylor coefficients of the ”normalized” Mather’s $\beta$–function are
invariants under $C^{\infty}$–conjugacies. In contrast, we prove that any two
elliptic billiard maps are $C^{0}$–conjugated near their respective
boundaries, and $C^{\infty}$–conjugated in the open cylinder, near the
boundary and away from a plain passing through the center of the underlying
ellipse. We also prove that if the billiard maps corresponding to two ellipses
are topologically conjugated then the two ellipses are similar.

MSC2020 numbers:  37C83, 37E40, 37J51

Keywords: Birkhoff billiard, integrability, conjugacy, Mather’s
${\beta}$–function, Marvizi-Melrose invariants.

##  1 Introduction

A billiard is a mathematical modeling of the dynamic of a confined massless
particle without friction and reflecting elastically on the boundary (without
friction): the particle moves along a straight line with constant speed till
it hits the boundary, then reflects off with reflection angle equals to the
angle of incidence and follows the reflected straight line. This seemingly
simple and fascinating Dynamical system was introduced by G.D. Birkhoff
[Bir20] in 1920. Since then, it has captured much of attention from both
Physicists and Mathematicians. We refer the reader to [Tab05] and references
therein for more details on Billiards. It is presently a very active and
popular subject with many challenging open problems.

The motivation of this work is the following question attributed to Victor
Guillemin:

###### Question 1

Let $f$ and $g$ be smooth Birkhoff billiard maps corresponding to two strictly
convex domains ${\Omega}_{f}$ and ${\Omega}_{g}$. Assume that f and g are
conjugate, i.e. there exists a homeomorphism h such that $f=h^{-1}\circ g\circ
h$. What can we say about the two domains ${\Omega}_{f}$ and ${\Omega}_{g}$?
Are they “similar”, that is, have they the same shape?

To our knowledge the only known answer to this question is in the case of
circular billiards. In fact, a billiard map in a disc $\mathscr{D}$ enjoys the
peculiar property of having the phase space completely foliated by
homotopically non-trivial invariant curves. It is an example of so–called
globally integrable billiards. In terms of the geometry of the billiard domain
$\mathscr{D}$, this reflects the existence of a smooth foliation by (smooth
and convex) caustics, i.e. (smooth and convex) curves with the property that
if a trajectory is tangent to one of them, then it will remain tangent after
each reflection. It is easy to check that in the circular billiard case this
family consists of concentric circles. See for instance [Tab05, Chapter 2].
The peculiarity of this property is justified by the following exciting result
due to Misha Bialy:

Theorem (M. Bialy[Bia93]) If the phase space of the billiard ball map is
foliated by continuous invariant curves which are not null–homotopic, then it
is a circular billiard.

In particular, the answer to Question 1 is affirmative whenever one of the
domains is a circle, since any foliation is invariant under conjugacy.

In general, billiards in ellipses are locally integrable: a neighborhood of
their boundary is foliated by caustics. Then, Birkhoff asked whether there are
other examples of (locally) integrable billiards. He actually conjectured:

######  Conjecture 1 (Birkhoff)

Amongst all convex billiards, only those in ellipses (circles being a distinct
special case) are integrable.

We are mainly concerned in this paper by invariant of Birkhoff billiards under
conjugacies when both of the domains are non–integrable and with the existence
of conjugacy (local and gloabal) for elliptic billiards tables (which
corresponds to the integrable case according to Birkhoff conjecture). before
summarizing the results of this work, let us introduce the notion of
”normalized” Mather’s $\beta$–function.

The Mather’s $\beta$–function (or Mather’s minimal average action) is a key
function in the celebrated Aubry–Mather Theory as it encodes several important
properties of the dynamics. The Mather’s $\beta$–function can be defined for
any exact area preserving twist map, not necessarily a billiard map. Roughly
speaking, its associates to any fixed rotation number (not only rational ones)
the minimal average action of orbits with that rotation number (whose
existence, inside a suitable interval, is ensured by the twist condition),
which we now define precisely.

For, let ${\Omega}$ be a bounded strictly convex domain in 2 with $C^{r}$
boundary ${\partial}{\Omega}$, with $r\geq 3$.333Observe that if ${\Omega}$ is
not convex, then the billiard map is not continuous; in this article we will
be interested only in strictly convex domains. Moreover, as pointed out by
Halpern [Hal77], if the boundary is not at least $C^{3}$, then the flow might
not be complete. The phase space $\mathscr{M}$ of the billiard map consists of
unit vectors $(x,v)$ whose foot points $x$ are on ${\partial}{\Omega}$ and
that have inward directions. The billiard ball map
$f:\mathscr{M}\to\mathscr{M}$ takes $(x,v)$ to $(x^{\prime},v^{\prime})$,
where $x^{\prime}$ represents the point where the trajectory starting at $x$
with velocity $v$ hits the boundary ${\partial}{\Omega}$ again, and
$v^{\prime}$ is the reflected velocity, according to the standard reflection
law: angle of incidence is equal to the angle of reflection. Assume that the
boundary ${\partial}{\Omega}$ is parametrized by arc–length $s$ and let
${\gamma}:[0,|{\partial}{\Omega}|]\to{}^{2}$ denote such a parametrization,
where $|{\partial}{\Omega}|$ denotes the length of ${\partial}{\Omega}$. Let
${\theta}$ be the angle between $v$ and the positive tangent to
${\partial}{\Omega}$ at $x$. Hence, $\mathscr{M}$ can be identified with the
annulus $\mathbb{A}\coloneqq[0,|{\partial}{\Omega}|]\times(0,\pi)$ and the
billiard map $f$ can be described as

$f\colon\mathbb{A}\ni(x,{\theta})\longmapsto(x^{\prime},{\theta}^{\prime})\in\mathbb{A}.$

Moreover, $f$ can be extended (smoothly) to
$\bar{\mathbb{A}}=[0,|{\partial}{\Omega}|]\times[0,\pi]$ by setting
$f(s,{\theta}_{0})=(s,{\theta}_{0})$, for all $s\in[0,|{\partial}{\Omega}|]$,
whenever ${\theta}_{0}\in\\{0,\pi\\}$. Denote the Euclidean distance between
two points of ${\partial}{\Omega}$ by

$d(s,s^{\prime})\coloneqq\|{\gamma}(s)-{\gamma}(s^{\prime})\|.$

Then, one checks easily

${\partial}_{s}d(s,s^{\prime})=-\cos{\theta},\qquad{\partial}_{s^{\prime}}d(s,s^{\prime})=\cos{\theta}^{\prime}.$
(1)

In particular, if we lift everything to the universal cover and introduce the
new coordinates $(x,y)=(s,\cos{\theta})\in\real\times[-1,1]$, then the
billiard map is a twist map with $d$ as generating function, and it preserves
the area form $dx\wedge dy$.444See [Tab05, Sib04] for details. Let
${\Gamma}\coloneqq\\{x_{n},y_{n}\\}_{n\geq 0}\subseteq\real\times[-1,1]$ be a
billiard trajectory.

###### Definition 2

If the limit $\lim_{n\to\pm\infty}x_{n}/n$ exists then it is called the
rotation number of the trajectory ${\Gamma}=\\{x_{n},y_{n}\\}_{n\geq 0}$ and
will be denoted by ${\rm\,rot\,}({\Gamma})$.

###### Remark 3

The trajectory ${\Gamma}\coloneqq\\{x_{n},y_{n}\\}_{n\geq 0}$ is called
$(p,q)$–periodic if $y_{q}=y_{0}$ and $x_{q}=x_{0}+p$. For such orbits, the
rotation number always exists and is equal to $p/q$. Furthermore, there always
exists at least two periodic orbits of rotation number $p/q$, for any rational
$p/q$ by a Theorem by Birkhoff [Bir13].

###### Definition 4

Given $p/q\in{\mathbb{Q}}$, denote by $\mathscr{P}er(p,q)$ the collection of
all the $(p,q)$–periodic trajectories. Thus, the Mather’s $\beta$–function at
$p/q$ is defined as

${\beta}(\frac{p}{q})\coloneqq-\max\\{{\,\rm length}({\Gamma})\ :\
{\Gamma}\in\mathscr{P}er(p,q)\\}.$

Then, ${\beta}$ is extended continuously to by setting

${\beta}({\omega})\coloneqq\lim_{n\to\infty}{\beta}(\frac{p_{n}}{q_{n}}),$

for any sequence $\\{{p_{n}}/{q_{n}}\\}_{n}$ of rationals converging to
${\omega}$.

###### Definition 5

For a given domain ${\Omega}$, we define the ”normalized” Mather’s
$\beta$–function as the function
${\lambda}^{-3}({\beta}+\ell\,{\mathtt{id}})$, where ${\beta}$ is the Mather’s
$\beta$–function, ${\lambda}$ is the Lazutkin perimeter, and $\ell$ the
perimeter of the domain ${\Omega}$.

The Taylor coefficients of the Mather’s $\beta$–function at $0$ are related to
the so-called Marvizi–Melrose invariants which was introduced by S. Marvizi
and R. Melrose in [MM82] where they proved the following. For, given
$q\in{\mathbb{N}}$, denote by

$\displaystyle L_{q}$ $\displaystyle\coloneqq\sup\\{{\,\rm length}({\Gamma})\
:\ {\Gamma}\in\mathscr{P}er(p,q),\ p\in{\mathbb{Z}}\\},$ $\displaystyle l_{q}$
$\displaystyle\coloneqq\inf\\{{\,\rm length}({\Gamma})\ :\
{\Gamma}\in\mathscr{P}er(p,q),\ p\in{\mathbb{Z}}\\},$

the supremum and the infimum of the perimeters of simple periodic billiard
trajectories with $q$ vertices, respectively. The:

###### Theorem 6

For any positive integer $k$ we have

$\lim_{q\to\infty}q^{k}(L_{q}-l_{q})=0.$

Moreover, $L_{q}$ has an asymptotic expansion as $q\to\infty$:

$L_{q}\sim\ell_{0}+\sum_{k=1}^{\infty}\frac{\ell_{k}}{q^{2k}},$

where $\ell_{0}$ is the length of the billiard table and $\ell_{k}$ are
constants, depending on the curvature of the table.

This collection $\\{\ell_{k}\\}_{k\geq 0}$ constitutes the Marvizi–Melrose
spectral invariants (also called Marvizi–Melrose invariants).

In the case of non–integrables domains, we prove that the Taylor coefficients
of the ”normalized” Mather’s $\beta$–function are invariant under smooth
conjugacies (see Theorem 7). As a conditional consequence, under Conjecture 8,
the sequence given by the ratios of the respective Marvizi–Melrose invariants
of such two $C^{\infty}$–conjugated domains is given by a geometric sequence
(see Corollary 9). In contrast, we prove that any two elliptic billiard maps
are $C^{0}$–conjugated near their respective boundaries, and
$C^{\infty}$–conjugated, in the open cylinder, near their respective
boundaries and off two lines. We also prove that the billiard maps
corresponding to two ellipses are topologically conjugated only if the two
ellipses are similar (i.e. they are the same up to a rescaling and an
isometry) or, equivalently, have the same eccentricity (see Theorem 11).555It
is an easy exercise to check that: two ellipses are similar iff they have the
same eccentricity. Finally, assuming the Birkhoff conjecture, we prove that
the answer to Question 1 is affirmative whenever one of the domains is
integrable (See Corollary 12).

##  2 Results

###### Theorem 7

Let ${\Omega}_{1},{\Omega}_{2}\subseteq{{}^{2}}$ be two strictly convex
domains with smooth boundaries. Let ${\beta}_{j}$ be the Mather’s
${\beta}$–function of ${\Omega}_{j}$, $\ell_{j}\coloneqq{\,\rm
length}({\partial}{\Omega}_{j})$, $\kappa_{j}(l)$ be the curvature of
${\partial}{\Omega}_{j}$ with arc–length parametrization $l$ and
${\lambda}_{j}\coloneqq\int_{0}^{\ell_{j}}\kappa_{j}(l)^{2/3}dl$.666Observe
that $\kappa_{j}>0$ by strictly convexity and, therefore, ${\lambda}_{j}>0$.
Assume that no strip in ${\Omega}_{2}\times[0,\pi)$ containing the boundary
${\Omega}_{2}\times\\{0\\}$ is foliated by invariant curves and

$f_{{\Omega}_{1}}\circ\tilde{h}=\tilde{h}\circ f_{{\Omega}_{2}},$ (2)

for some homeomorphism $\tilde{h}\in C^{\infty}({\mathbb{T}}\times\real)$.
Then, ${\lambda}_{1}^{-3}({\beta}_{1}({\omega})+\ell_{1}{\omega})$ and
${\lambda}_{2}^{-3}({\beta}_{2}({\omega})+\ell_{2}{\omega})$ have the same
Taylor’s expansion at ${\omega}=0$. In particular,
${\lambda}_{1}^{-3}\mathscr{L}_{1}({\omega})$ and
${\lambda}_{2}^{-3}\mathscr{L}_{2}({\omega})$ have the same Taylor’s expansion
at ${\omega}=0$, where $\mathscr{L}_{j}({\omega})$ is the Lazutkin’s parameter
of ${\Omega}_{j}$ corresponding to the convex caustic of rotation number
${\omega}$ whenever it exists.

Denote by by $\\{\mathcal{I}^{j}_{2n+1}\\}_{n}$ the Marvizi–Melrose invariants
of ${\Omega}_{j}$. We conjecture the following relation between the Taylor’s
coefficients of the Mather’s ${\beta}$–function and the Marvizi–Melrose
invariants.777The Marvizi–Melrose invariants are indeed known to be
algebraically equivalent to the Taylor’s coefficients of the Mather’s
${\beta}$–function but, to our best knowledge, no explicit algebraic formula
relating them in known.

###### Conjecture 8

For any $n\in{\mathbb{N}}$,

$\frac{d^{2n+1}{\beta}}{d{\omega}^{2n+1}}(0)=\mathcal{I}_{3}^{n+2}\cdot\sum_{\begin{subarray}{c}{\sigma}_{1},\cdots,{\sigma}_{n}\in{\mathbb{N}}\cup\\{0\\}\\\
{\sigma}_{1}+\cdots+{\sigma}_{n}=n-1\\\
3{\sigma}_{1}+\cdots+(2n+1){\sigma}_{n}=5(n-1)\end{subarray}}r_{n}({\sigma}_{1},\cdots,{\sigma}_{n})\cdot\mathcal{I}_{3}^{{\sigma}_{1}}\cdots\mathcal{I}_{2n+1}^{{\sigma}_{n}}\;,$
(3)

where each $r_{n}({\sigma}_{1},\cdots,{\sigma}_{n})\in{\mathbb{Q}}$.

Under Conjecture 8, Theorem 7 yields easily (e.g. by induction) the following.

###### Corollary 9

The sequence $\\{\mathcal{I}^{1}_{2n+1}/\mathcal{I}^{2}_{2n+1}\\}_{n}$ is a
geometric sequence:888Observe that ${\lambda}_{j}=\mathcal{I}^{j}_{3}$.

$\mathcal{I}^{1}_{2n+1}=({\lambda}_{2}/{\lambda}_{1})^{2n-3}\;\mathcal{I}^{2}_{2n+1},\qquad\forall\;n\geq
1.$ (4)

###### Remark 10

Using [Sor15], one checks easily (3) and (4) up to $n=4$.

In contrast to Theorem 7, the following holds for elliptic Billiards and,
hence, according to Birkhoff Conjecture 8, for convex integrable billiards.

###### Theorem 11

Let ${\Omega}_{j}$ be an ellipse with eccentricity $\mathbf{e}_{j}$ and
$f_{j}$ be the associated billiard map ($j=1,2$). Set999We use $-p$ to denote
the symmetric of the point $p$ w.r.t the origin.

${\Omega}_{j}^{*}\coloneqq{\Omega}_{j}\setminus\\{p_{j}^{*},-p_{j}^{*}\\},\quad\mbox{where}\quad
p_{j}^{*}\coloneqq(a_{j}\cos\arctan^{-1}(1-\mathbf{e}_{j}^{2}),-b_{j}\sin\arctan^{-1}(1-\mathbf{e}_{j}^{2})).$

Then

$(i)$ The billiard maps $f_{1}$ and $f_{2}$ are $C^{0}$–conjugated near their
respective boundaries i.e., there exists ${\theta}_{j}^{*}>0$ such that their
restrictions $f_{1}|_{{\Omega}_{1}\times(0,{\theta}_{1}^{*})}$ and
$f_{2}|_{{\Omega}_{2}\times(0,{\theta}_{2}^{*})}$ are $C^{0}$–conjugated.
Moreover, the restrictions
$f_{1}|_{{\Omega}_{1}^{*}\times(0,{\theta}_{1}^{*})}$ and
$f_{2}|_{{\Omega}_{2}^{*}\times(0,{\theta}_{2}^{*})}$ are
$C^{\infty}$–conjugated.

$(ii)$ If the billiard maps $f_{1}$ and $f_{2}$ are (globally)
$C^{0}$–conjugated then ${\Omega}_{1}$ and ${\Omega}_{2}$ are similar (or,
equivalently, $\mathbf{e}_{1}=\mathbf{e}_{2}$).

As a consequence, assuming the Birkhoff conjecture, we have:

###### Corollary 12

Let $f$ and $g$ be smooth Birkhoff billiard maps corresponding to two strictly
convex domains ${\Omega}_{f}$ and ${\Omega}_{g}$. Assume that $f$ and $g$ are
topologically conjugated and that one of the two domains is integrable. Then,
${\Omega}_{f}$ and ${\Omega}_{g}$ are similar.

##  3 Proofs

###  3.1 Proof of Theorem 7

Given $j=1,2$, pick a symplectic change of coordinates $\phi_{j}$ with
associated Cantor set $\mathscr{C}_{j}$ corresponding to $f_{{\Omega}_{j}}$
given by Theorem A.2; set

$f_{j}\coloneqq\phi_{j}^{-1}\circ
f_{{\Omega}_{j}}\circ\phi_{j},\quad\mbox{and}\quad
h=(h_{1},h_{2})\coloneqq\phi_{1}^{-1}\circ\tilde{h}\circ\phi_{2}.$

Let ${\beta}_{j}\coloneqq{\beta}_{{\Omega}_{j}}$, the minimal average action
(or Mather’s ${\beta}$–function) of $f_{{\Omega}_{j}}$. Observe that (2) reads

$f_{1}\circ h=h\circ f_{2},$ (5)

which, recalling101010We refer the reader to [Sor15] for explicit formula of
the coefficients of the Taylor’s expansion of the ${\alpha}$ and ${\beta}$
functions at the boundary.
${\alpha}_{1}^{\prime}(-\ell_{1})={\alpha}_{1}^{\prime}(-\ell_{1})=0$, implies
for ${\theta}=-\ell_{2}$,

${\alpha}_{1}^{\prime}(h_{2}(s,-\ell_{2}))=0={\alpha}_{1}^{\prime}(-\ell_{1}),\qquad\forall\;s\in{\mathbb{T}}.$

Thus,

$h_{2}(s,-\ell_{2})=-\ell_{1},\qquad\forall\;s\in{\mathbb{T}}.$ (6)

Given $j=1,2$ and ${\omega}\in{\alpha}_{j}^{\prime}(\mathscr{C}_{j})$, let
denote by
${\Gamma}^{j}_{\omega}\coloneqq{\mathbb{T}}\times\\{{\beta}_{j}^{\prime}({\omega})\\}$,
the KAM torus of frequency ${\omega}$ of $f_{j}$ and by
$\mathscr{D}^{j}_{\omega}$, the domain of the cylinder
${\mathbb{T}}\times\real$ bounded by the closed curves
$\\{{\theta}={\beta}_{j}({\omega})\\}$ and $\\{{\theta}=-\ell_{j}\\}$.

Combining (5) and (A.2), we get that, for any
${\omega}\in{\alpha}_{2}^{\prime}(\mathscr{C}_{2})\setminus\\{0\\}$,
$h({\Gamma}^{2}_{\omega})$ is an invariant curve of Diophantine frequency
${\omega}$ for $f_{1}$. Thus, by Lemma B.1, $h({\Gamma}^{2}_{\omega})$ is
Lagrangian and, hence, is a KAM curve of frequency ${\omega}$ of $f_{1}$.
Therefore, by uniqueness, for any
${\omega}\in{\alpha}_{2}^{\prime}(\mathscr{C}_{2})\setminus\\{0\\}$,

$h({\Gamma}^{2}_{\omega})={\Gamma}^{1}_{\omega},$ (7)

and, by continuity, for any
$(s,{\theta})\in{\mathbb{T}}\times\mathscr{C}_{2}$,

${\alpha}_{1}^{\prime}(h_{2}(s,{\theta}))={\alpha}_{2}^{\prime}({\theta}).$
(8)

Thus, as $h$ is a homeomorphism of the cylinder ${\mathbb{T}}\times\real$, we
infer

$h(\mathscr{D}^{2}_{\omega})=\mathscr{D}^{1}_{\omega}\qquad\mbox{and}\qquad
h({\partial}\mathscr{D}^{2}_{\omega})={\partial}\mathscr{D}^{1}_{\omega}.$ (9)

Set ${\,\rm J}h\coloneqq\det{\partial}_{(s,{\theta})}h$.111111 Here and
thereafter, we denote by ${\partial}_{(s,{\theta})}h$ the Jacobian matrix of
$h$. We claim that

${\partial}_{\theta}^{k}{\,\rm
J}h(s,-\ell_{2})=0,\qquad\forall\;s\in{\mathbb{T}},\ \forall\;k\geq 1.$ (10)

Fix $s\in{\mathbb{T}}$. Indeed, we are going to construct inductively, for a
given $k\geq 1$, a sequence $\\{{\omega}_{n}^{k}\\}_{n}\searrow 0$ such that

${\partial}_{\theta}^{k}{\,\rm
J}h(s,{\beta}_{2}^{\prime}({\omega}_{n}^{k}))=0,\qquad\forall\;n\geq 1,$ (11)

which implies (10) by continuity. For, pick a strictly decreasing sequence
$\\{{\omega}_{n}^{0}\\}_{n}\subseteq{\alpha}_{2}^{\prime}(\mathscr{C}_{2})\setminus\\{0\\}$
converging to $0$.121212Such sequence exists as
${\alpha}_{2}^{\prime}(\mathscr{C}_{2})$ is a perfect set. As no strip near
the boundary of ${\Omega}_{2}$ is foliated by invariant curves, then, by
Mather connecting Theorem (see [Mat90, Kal03]) any two successive invariant
curves of $f_{2}$ are (asymptotically) connected by some orbit, the invariant
curves being naturally ordered via their respective rotation numbers. Let then
$\mathscr{M}_{2}$ be the closure of the union of all invariant curves and
Mather’s connecting orbits of $f_{2}$; in particular, $\mathscr{M}_{2}$
contains the family of KAM tori of $f_{2}$. As $f_{1}$ and $f_{2}$ are
symplectic, we have
$\det{\partial}_{(s,{\theta})}f_{1}=\det{\partial}_{(s,{\theta})}f_{2}\equiv
1$ and, therefore, differentiating (5), we obtain, for any
$(s,{\theta})\in{{}^{2}}$, ${\,\rm J}h(f_{2}(s,{\theta}))={\,\rm
J}h(s,{\theta})$ i.e. ${\,\rm J}h$ is constant along orbits of $f_{2}$.
Consequently, ${\,\rm J}h$ is constant on $\mathscr{M}_{2}$ by continuity.
Then, for any $n\in{\mathbb{N}}$, as ${\,\rm
J}h(s,{\beta}_{2}^{\prime}({{\omega}}_{n}^{0}))={\,\rm
J}h(s,{\beta}_{2}^{\prime}({{\omega}}_{n+1}^{0}))$, by the Mean Value Theorem,
there exists ${\omega}_{n+1}^{0}<{{\omega}}_{n}^{1}<{\omega}_{n}^{0}$ such
that ${\partial}_{\theta}{\,\rm
J}h(s,{\beta}_{2}^{\prime}({{\omega}}_{n}^{1}))=0.$ Observe that, by
construction, $\\{{\omega}_{n}^{1}\\}_{n}\searrow 0$, which proves (11) for
$k=1$. Moreover, we recover the setting we started with, with now
$\\{{\omega}_{n}^{1}\\}_{n}$ playing the role of $\\{{\omega}_{n}^{0}\\}_{n}$.
Iterating the argument, we then obtain (11) and, in particular, (10). Hence,
for ${\omega}$ near $0$,131313Recall that
${\beta}_{2}^{\prime}({\omega})+\ell_{2}=O({\omega}^{2})$; see [Sor15].

${\,\rm
J}h(s,{\beta}_{2}^{\prime}({\omega}))\stackrel{{\scriptstyle{\rm(\ref{dthhOinfo})}}}{{=}}{\,\rm
J}h(s,-\ell_{2})+O(({\beta}_{2}^{\prime}({\omega})+\ell_{2})^{k})={\,\rm
J}h(s,-\ell_{2})+O({\omega}^{2k}),\qquad\forall\;k\geq 1.$ (12)

Thus, for any ${\omega}\in{\alpha}_{1}^{\prime}(\mathscr{C}_{1})$ and $k\geq
1$,

$\displaystyle{\beta}_{1}^{\prime}({\omega})+\ell_{1}=\int_{{\partial}\mathscr{D}_{\omega}^{1}}{\theta}ds$
$\displaystyle=\int_{\mathscr{D}_{\omega}^{1}}d{\theta}ds\qquad\mbox{(by
Stokes--Whitney's Theorem)}$
$\displaystyle=\int_{h(\mathscr{D}_{\omega}^{1})}|{\,\rm
J}h|d{\theta}ds\stackrel{{\scriptstyle{\rm(\ref{boudry})}}}{{=}}\int_{\mathscr{D}_{\omega}^{2}}|{\,\rm
J}h|d{\theta}ds\quad\mbox{(by making
$(s^{\prime},{\theta}^{\prime})=h(s,{\theta})$)}$
$\displaystyle=\int_{{\partial}\mathscr{D}_{\omega}^{2}}|{\,\rm
J}h|{\theta}ds\qquad\qquad\qquad\qquad\quad\mbox{(by Stokes--Whitney's
Theorem)}$ $\displaystyle=\int_{{\mathbb{T}}}(|{\,\rm
J}h(s,{\beta}^{\prime}_{2}({\omega}))|{\beta}^{\prime}_{2}({\omega})+|{\,\rm
J}h(s,-\ell_{2})|\ell_{2})\;ds$
$\displaystyle\stackrel{{\scriptstyle{\rm(\ref{oDInf})}}}{{=}}\int_{{\mathbb{T}}}(({\beta}_{2}^{\prime}({\omega})+\ell_{2})|{\,\rm
J}h(s,-\ell_{2})|+O({\omega}^{2k}))\;ds$
$\displaystyle=({\beta}_{2}^{\prime}({\omega})+\ell_{2})\int_{{\mathbb{T}}}|{\,\rm
J}h(s,-\ell_{2})|ds+O({\omega}^{2k})$ (13)

We now compute $\int_{{\mathbb{T}}}|{\,\rm J}h(s,-\ell_{2})|ds$. In fact, by
Lemma 13, we have ${\,\rm J}h(s,{\theta})\equiv h_{2}^{\prime}({\theta})$.
But141414See [Sor15] .

$({\theta}+\ell_{j})^{-1/2}{\alpha}_{j}^{\prime}({\theta})=2\sqrt{2}{\lambda}_{j}^{-3/2}+O({\theta}+\ell_{j}).$

Thus, rewriting the second part of (15) as

$\left(\frac{h_{2}({\theta})+\ell_{1}}{{\theta}+\ell_{2}}\right)^{1/2}(h_{2}({\theta})+\ell_{1})^{-1/2}{\alpha}_{1}^{\prime}(h_{2}({\theta}))=({\theta}+\ell_{2})^{-1/2}{\alpha}_{2}^{\prime}({\theta}),$

and using (6), it follows

${\,\rm J}h(s,-\ell_{2})\equiv
h_{2}^{\prime}(-\ell_{2})=\left(\frac{{\lambda}_{1}}{{\lambda}_{2}}\right)^{3}.$

Hence, for any $k\geq 1$,

${\lambda}_{1}^{-3}({\beta}_{1}^{\prime}({\omega})+\ell_{1})={\lambda}_{2}^{-3}({\beta}_{2}^{\prime}({\omega})+\ell_{2})+O({\omega}^{2k}),$

and, therefore,

${\lambda}_{1}^{-3}({\beta}_{1}({\omega})+\ell_{1}{\omega})={\lambda}_{2}^{-3}({\beta}_{2}({\omega})+\ell_{2}{\omega})+O({\omega}^{2k+1}),$
(14)

i.e. ${\lambda}_{1}^{-3}({\beta}_{1}({\omega})+\ell_{1}{\omega})$ and
${\lambda}_{2}^{-3}({\beta}_{2}({\omega})+\ell_{2}{\omega})$ have the same
Taylor expansion. Now, for ${\omega}$ small enough, if a convex caustic of
rotation number ${\omega}$ of ${\Omega}_{j}$ exists, then it Lazutkin’s
parameter is given by
$\mathscr{L}_{j}({\omega})={\omega}{\beta}_{j}^{\prime}({\omega})-{\beta}_{j}({\omega})$
and, therefore,
${\lambda}_{1}^{-3}\mathscr{L}_{1}({\omega})={\lambda}_{2}^{-3}\mathscr{L}_{2}({\omega})+O({\omega}^{2k+1})$.

###### Lemma 13

Under the notations and assumptions above, there exists151515Actually, $g_{1}$
is as smooth as $h$. $g_{1}\in C^{\infty}(\real,\real)$ such that, for any
$(s,{\theta})\in{\mathbb{T}}\times\mathscr{C}_{2}$, with ${\theta}$
sufficiently close to $-\ell_{2}$,

$h(s,{\theta})=(g_{1}({\theta})+s,h_{2}({\theta})),\qquad{\alpha}_{1}^{\prime}(h_{2}({\theta}))={\alpha}_{2}^{\prime}({\theta}).$
(15)

Proof Let
$(s,{\theta})\in{\mathbb{T}}\times\mathscr{C}_{2}\setminus\\{-\ell_{2}\\}$.
Then, we have

$\displaystyle\big{(}h_{1}(s,{\theta})+{\alpha}_{1}^{\prime}(h_{2}(s,{\theta})\big{)},h_{2}(s,{\theta}))\stackrel{{\scriptstyle{\rm(\ref{conjBill})}}}{{=}}f_{1}\circ
h(s,{\theta})\stackrel{{\scriptstyle{\rm(\ref{conjf1f2})}}}{{=}}h\circ
f_{2}(s,{\theta})\stackrel{{\scriptstyle{\rm(\ref{conjBill})}}}{{=}}$
$\displaystyle\qquad\stackrel{{\scriptstyle{\rm(\ref{conjBill})}}}{{=}}\big{(}h_{1}(s+{\alpha}_{2}^{\prime}({\theta}),{\theta}),h_{2}(s+{\alpha}_{2}^{\prime}({\theta}),{\theta})\big{)}.$
(16)

Hence, $h_{2}(s,{\theta})=h_{2}(s+{\alpha}_{2}^{\prime}({\theta}),{\theta})$
which implies $h_{2}$ is independent of $s$ i.e.
$h_{2}(s,{\theta})=h_{2}({\theta})$, as
${\alpha}_{2}^{\prime}({\theta})\not\in{\mathbb{Q}}$. We have also

$h_{1}(s,{\theta})+{\alpha}_{1}^{\prime}(h_{2}({\theta}))\stackrel{{\scriptstyle{\rm(\ref{deteH})}}}{{=}}h_{1}(s+{\alpha}_{2}^{\prime}({\theta}),{\theta}),$
(17)

which differentiated w.r.t $s$ yields
${\partial}_{s}h_{1}(s,{\theta})={\partial}_{s}h_{1}(s+{\alpha}_{2}^{\prime}({\theta}),{\theta})$.
Thus, as ${\alpha}_{2}^{\prime}({\theta})\not\in{\mathbb{Q}}$,
${\partial}_{s}h_{1}$ is also independent of $s$:
${\partial}_{s}h_{1}(s,{\theta})=g_{2}({\theta})$, so that
$h_{1}(s,{\theta})=g_{1}({\theta})+g_{2}({\theta})\cdot s$, for some
$g_{1},g_{2}\in C^{\infty}(\real,\real)$. Now, plugging the expression of
$h_{1}$ found into (17) and using (A.1), we get

$h(s,{\theta})=(g_{1}({\theta})+g_{2}({\theta})\cdot
s,h_{2}({\theta})),\qquad{\alpha}_{1}^{\prime}(h_{2}({\theta}))=g_{2}({\theta}){\alpha}_{2}^{\prime}({\theta}).$
(18)

Now, combining (8), (18) with the fact that161616See [Sor15].
${\alpha}_{2}^{\prime}({\theta})\not=0$ for ${\theta}+\ell_{2}>0$ small enough
yields $g_{2}({\theta})=1,$ for ${\theta}$ near $-\ell_{2}$. Thus, as
$\mathscr{C}_{2}\setminus\\{-\ell_{2}\\}$ accumulates at $-\ell_{2}$, by
continuity we obtain $g_{2}\equiv 1$, and (15) is proven.

###### Remark 14

(i) Observe that, according to the proof, the conclusion in Lemma 13 still
holds if $\tilde{h}$ is assumed merely to be
$C^{2}({\mathbb{T}}\times\real,\real)$.
(ii) Observe also that, if ${\Omega}_{1}$ and ${\Omega}_{2}$ are both ellipses
and $\tilde{h}$ is $C^{2}$ conjugacy near their respective boundaries of the
associated billiard maps $f_{{\Omega}_{1}}$ and $f_{{\Omega}_{2}}$, then (7)
and (8), and, therefore, (15) hold is some neighborhood
$[-\ell_{2},-\ell_{2}+{\varepsilon}_{2})$ of $-\ell_{2}$. Indeed, by the same
argument, one gets that (7), (8) and (15) hold on
$[-\ell_{2},-\ell_{2}+{\varepsilon}_{2})\setminus{\mathbb{Q}}$ and, therefore,
by continuity they hold on $[-\ell_{2},-\ell_{2}+{\varepsilon}_{2})$.

###  3.2 Proof of Theorem 11

####  3.2.1 Proof of $(i)$ in Theorem 11

We shall exhibit a $C^{\infty}$–conjugacy. Given $j=1,2$ let $0<b_{j}\leq
a_{j}$ such that

${\Omega}_{j}\coloneqq\left\\{(x,y)\in{{}^{2}}\ :\
\frac{x^{2}}{a_{j}^{2}}+\frac{y^{2}}{b_{j}^{2}}=1\right\\},$

so that $\mathbf{e}_{j}\coloneqq\sqrt{1-(b_{j}/a_{j})^{2}}$. Let
$\mathsf{c}_{j}\coloneqq\sqrt{a_{j}^{2}-b_{j}^{2}}$ its semi–focal distance,
and $\mathcal{F}_{j}^{\pm}\coloneqq(\pm\mathsf{c}_{j},0)$ its two foci.
Consider the family of confocal elliptic caustics

${\Omega}_{j,{\lambda}}\coloneqq\left\\{(x,y)\in{{}^{2}}\ :\
\frac{x^{2}}{a_{j}^{2}-{\lambda}^{2}}+\frac{y^{2}}{b_{j}^{2}-{\lambda}^{2}}=1\right\\}\;,\qquad{\lambda}\in[0,b).$

We denote the (incomplete) elliptic integral of the first kind by

$F\colon[0,2\pi)\times[0,1)\ni({\varphi},k)\longmapsto\int_{0}^{\varphi}\frac{d\phi}{\sqrt{1-k^{2}\sin^{2}\phi}},$

and the complete elliptic integral of the first kind by $K(k)\coloneqq
F(\pi/2,k)$.

Consider the following auxiliary functions

$\displaystyle
k_{j}\colon[0,b_{j})\ni{\lambda}\longmapsto\sqrt{\frac{a_{j}^{2}-b_{j}^{2}}{a_{j}^{2}-{\lambda}^{2}}}\in[0,1)\;,\qquad{\omega}_{j}\colon[0,b_{j})\ni{\lambda}\longmapsto\frac{F(\arcsin({\lambda}/b_{j}),k_{j}({\lambda}))}{2K(k_{j}({\lambda}))}\in[0,\frac{1}{{2}})\;,$
$\displaystyle
h_{j,1}\colon[0,2\pi)\times[0,\pi)\ni(\phi,{\theta})\longmapsto(a_{j}\cos\phi,b_{j}\sin\phi,{\theta})\in{\Omega}_{j}\times[0,\pi),$
$\displaystyle
h_{j,2}\colon[0,2\pi)\times[0,\pi)\ni(\phi,{\theta})\longmapsto({\lambda},t)\coloneqq\bigg{(}{g_{j}(\phi,{\theta})},F\big{(}\phi,k_{j}({g_{j}(\phi,{\theta})}\big{)}\bigg{)}\in[0,b_{j})\times[0,4K(k_{j}(b_{j})))$
$\displaystyle
h_{j,3}\colon[0,1)\times[0,1)\ni({\lambda},t)\longmapsto\big{(}b_{j}{\lambda},4K(k_{j}(b_{j}{\lambda}))t\big{)}\in[0,b_{j})\times[0,4K(k_{j}(b_{j})))\;,$

where

$\displaystyle g_{j}(\phi,{\theta})\coloneqq$
$\displaystyle\frac{|a_{j}^{2}\sin\phi+b_{j}^{2}\cos\phi|\sin{\theta}}{\sqrt{a_{j}^{2}\sin^{2}\phi+b_{j}^{2}\cos^{2}\phi}}\;.$
(19)

It is obvious that each of the maps $k_{j}$, ${\omega}_{j}$, $h_{j,1}$ and
$h_{j,3}$ are $C^{\infty}$–diffeomorphism (onto their respective codomains).
The map $h_{j,2}$ is a homeomorphism (onto its codomain) and its restriction
$h_{j,2}^{*}\colon\mathscr{T}^{*}\times(0,\pi)\overset{onto}{\longrightarrow}(0,b_{j})\times[0,4K(k_{j}(b_{j})))$
is a $C^{\infty}$–diffeomorphism, where
$\mathscr{T}^{*}\coloneqq[0,2\pi)\setminus\\{-\arctan^{-1}(b_{j}^{2}/a_{j}^{2}),\pi-\arctan^{-1}(b_{j}^{2}/a_{j}^{2})\\}$.171717Indeed,
the set of singularities of $h_{j,2}$ is precisely
$g_{j}^{-1}(\\{0\\})=\\{(\phi,{\theta})\in[0,2\pi)\times[0,\pi):\
\tan\phi=-b_{j}^{2}/a_{j}^{2}\mbox{ or }{\theta}=0\\}$ In particular, the
following map is a $C^{\infty}$–diffeomorphism as well:

$h_{j,4}\colon[0,1)\times[0,1)\ni({\lambda},t)\longmapsto({\omega}_{1}^{-1}\circ{\omega}_{2}({\lambda}),t)\in[0,1)\times[0,1).$

Given $s=(x,y)\in{\Omega}_{j}$, let181818Here and thereafter, $[P,Q]$ denotes
the closed segment joining $P$ and $Q$ i.e.
$[P,Q]\coloneqq\\{{\alpha}P+(1-{\alpha})Q\;:\;0\leq{\alpha}\leq 1\\}$.

$\displaystyle{\theta}_{j}(s)$
$\displaystyle\coloneqq\sup\left\\{{\varepsilon}\in(0,\pi)\ :\
\bigcup_{n=0}^{\infty}[f_{j}^{n}(s,{\theta}),\;f_{j}^{n+1}(s,{\theta})]\bigcap[\mathcal{F}_{j}^{-},\;\mathcal{F}_{j}^{+}]=\emptyset\;,\
\forall\;{\theta}\in[0,{\varepsilon}]\right\\}\;.$

Then, one checks easily that

${\theta}_{j}(x,0)=\frac{\pi}{2}\qquad\mbox{and}\qquad{\theta}_{j}(x,y)=\arctan\frac{b_{j}^{2}}{\mathsf{c}_{j}|y|},\quad\mbox{for}\
y\not=0.$

In particular,

${\theta}_{j}^{*}\coloneqq\inf_{s\in{\Omega}_{j}}{\theta}_{j}(s)=\arctan\frac{b_{j}^{2}}{\mathsf{c}_{j}}>0.$

Let ${\theta}_{3}^{*}\in[0,\pi)$ such that
$h^{-1}({\Omega}_{1}\times[0,{\theta}_{1}^{*}))={\Omega}_{2}\times[0,{\theta}_{3}^{*})$
and define

${\theta}^{*}\coloneqq\min\\{{\theta}_{2}^{*},{\theta}_{3}^{*}\\}.$

Furthermore, consider

$h_{j}\coloneqq h_{j,1}\circ h_{j,2}^{-1}\circ h_{j,3}\circ
h_{j,4}:[0,1)\times[0,1)\to{\Omega}_{j}\times[0,2\pi)$

and

$h\coloneqq h_{1}\circ
h_{2}^{-1}\colon{\Omega}_{2}\times[0,\pi)\to{\Omega}_{1}\times[0,\pi).$

Then, $h$ is a homeomorphism and its restriction

$h^{*}\colon{\Omega}_{2}^{*}\times(0,\pi)\to{\Omega}_{1}^{*}\times(0,\pi).$

is a $C^{\infty}$–diffeomorphisms. Now, we claim that

$f_{1}\circ h=h\circ
f_{2}\qquad\mbox{on}\qquad{\Omega}_{2}\times[0,{\theta}^{*}),$ (20)

which would prove $(i)$ in Theorem 11. The prove of (20) follows easily a
result proven in [CF88] (see also [KS18, Proposition 16]), which can be
reformulated as follows. The main difference with the aforementioned result is
that we provide the explicit value of $({\lambda},t)$, the “action–angle”
coordinates (see Appendix C for the computation of ${\lambda}$).

###### Lemma 15

Let $(s,{\theta})\in{\Omega}_{j}\times[0,\pi)$. Set $(\phi,{\theta})\coloneqq
h_{j,1}^{-1}(s,{\theta}),\ ({\lambda},t)\coloneqq h_{j,2}(\phi,{\theta})$ and
${\delta}_{j}({\lambda})\coloneqq
2F(\arcsin({\lambda}/b_{j}),k_{j}({\lambda}))=4K(k_{j}({\lambda})){\omega}_{j}({\lambda}).$
Then, the segment joining $(s,{\theta})$ and $f_{j}(s,{\theta})$ is tangent to
${\Omega}_{j,{\lambda}}$ and

$f_{j}\circ h_{j,1}\circ
h_{j,2}^{-1}({\lambda},t)=f_{j}(s,{\theta})=h_{j,1}\circ
h_{j,2}^{-1}({\lambda},t+{\delta}_{j}({\lambda})).$ (21)

####  3.2.2 Proof of $(ii)$ in Theorem 11

We adopt the notations in §3.2.1. Given $j=1,2$ and two coprime integers
$0<m<n/2$, consider the function

$g_{j}\colon(-\mathsf{c}_{j}^{2},0)\ni{\xi}\longmapsto
F\left(\arcsin\sqrt{\frac{b_{j}^{2}}{b_{j}^{2}-{\xi}}},\sqrt{1+\frac{{\xi}}{c_{j}^{2}}}\right)-\frac{2m}{n}F\left(\frac{\pi}{2},\sqrt{1+\frac{{\xi}}{c_{j}^{2}}}\right).$

Then191919See [Sie97, pg. 4569, Eq. (22)]., $f_{j}$ has a periodic orbit of
rotation $m/n$ whose caustic is a hyperbola iff the equation

$g_{j}({\xi})=0,\qquad{\xi}\in(-\mathsf{c}_{j}^{2},0)$ (22)

has a solution. But, we have

$g_{j}(-\mathsf{c}_{j}^{2})=\arcsin\frac{b}{a}-\frac{m}{n}\pi,\qquad\mbox{and}\qquad\lim_{{\xi}\to
0^{-}}g_{j}({\xi})=+\infty.$

Hence, if

$\frac{m}{n}\geq\frac{1}{\pi}\arcsin\frac{b_{j}}{a_{j}},$ (23)

then (22) has a solution. In contrast, if

$\frac{m}{n}<\frac{1}{\pi}\arcsin\frac{b_{j}}{a_{j}},$ (24)

then (22) has no solution. Indeed, assuming (24), we have, for any
${\xi}\in(-\mathsf{c}_{j}^{2},0)$,

$g_{j}({\xi})>u_{j}({\xi})>0,$

where

$u_{j}({\xi})\coloneqq
F\left(\arcsin\sqrt{\frac{b_{j}^{2}}{b_{j}^{2}-{\xi}}},\sqrt{1+\frac{{\xi}}{c_{j}^{2}}}\right)-\frac{2}{\pi}\arcsin(\frac{b_{j}}{a_{j}})F\left(\frac{\pi}{2},\sqrt{1+\frac{{\xi}}{c_{j}^{2}}}\right).$

and, using Wolfram Matematica[Inc], one gets

$\min\left\\{u_{j}({\xi}):\ a_{j}>b_{j}>0,\
-a_{j}^{2}+b_{j}^{2}<{\xi}<0\right\\}=5.65246\cdot 10^{-9}.$

Now, assume $f_{1}$ and $f_{2}$ are (globally) $C^{0}$–conjugated. By
contradiction, assume $\mathsf{e}_{1}>\mathsf{e}_{2}$ i.e.
${b_{1}}/{a_{1}}<{b_{2}}/{a_{2}}$. Then, for any two coprime integers
$0<m<n/2$ such that,

$\frac{1}{\pi}\arcsin\frac{b_{1}}{a_{1}}\leq\frac{m}{n}<\frac{1}{\pi}\arcsin\frac{b_{2}}{a_{2}}\,,$

$f_{1}$ admits a periodic orbit of rotation $m/n$ whose caustic is a hyperbola
while $f_{2}$ does not admit such a periodic orbit, contradiction. Thus,
$\mathsf{e}_{1}=\mathsf{e}_{2}.$

## Appendix

##  A Auxiliary facts about billiard dynamics

######  Theorem A.1 ([MF94, Mat90, Sib04])

Let $f$ be a monotone twist map. Then:

  * $(i)$

${\beta}$ is strictly convex. In particular, it is continuous and admits a
Legendre–Fenchel transform

${\alpha}(c)=\sup_{{\omega}\in\real}{\omega}\cdot c-{\beta}({\omega}),$

the so–called Mather’s ${\alpha}$ function. Moreover, ${\alpha}$ is convex and
we have

${\alpha}^{\prime}({\beta}^{\prime}({\omega}))={\omega},$ (A.1)

at each point of differentiability ${\omega}$ of ${\beta}$.

  * $(ii)$

${\beta}$ is differentiable at any irrational.

  * $(iii)$

${\beta}$, and hence ${\alpha}$, is an invariant under any symplectic change
of coordinates.

######  Theorem A.2 ([KP90, Sib04])

Let ${\Omega}\subseteq{{}^{2}}$ be a smooth strictly convex closed curve of
length $\ell$ and $f$ its corresponding billiard map. Then, there exists
${\varepsilon}_{0}>0$ small, a cantor set
$\mathscr{C}\subseteq[-\ell,-\ell+{\varepsilon}_{0})$, with
$-\ell\in\mathscr{C}$, a smooth (exact) symplectomorphism
$\phi\colon{\mathbb{T}}\times[-\ell,-\ell+{\varepsilon}_{0})\ni(x_{0},{\theta}_{0})\longmapsto(x_{1},{\theta}_{1})\in\phi({\mathbb{T}}\times[-\ell,-\ell+{\varepsilon}_{0}))\subseteq{\mathbb{T}}\times\real$
such that

$\phi^{-1}\circ
f\circ\phi(x,{\theta})=(x+{\alpha}^{\prime}({\theta}),{\theta}),\qquad\forall\;(x,{\theta})\in{\mathbb{T}}\times\mathscr{C}.$
(A.2)

In particular, the KAM curve of rotation number
${\omega}\in{\alpha}^{\prime}(\mathscr{C})$ of $\phi^{-1}\circ f\circ\phi$ is
given by the graph of the constant function ${\mathbb{T}}\ni
x\longmapsto{\beta}^{\prime}({\omega})$.

##  B Isotropicity Lemma

######  Lemma B.1 ([Her89, BHS09])

Let $f\colon\mathcal{M}\coloneqq{}^{m}\times{\mathbb{T}}^{m}\to
f(\mathcal{M})\subseteq\mathcal{M}$ be a symplectic202020w.r.t the canonical
symplectic form. diffeomorphism of class $C^{1}$and
$\phi:{\mathbb{T}}^{m}\to\mathcal{M}$, a $C^{1}$–embedding. Assume that

$f\circ\phi\circ g(x)=\phi\circ
g(x+t{\omega}),\qquad\forall\;x\in{\mathbb{T}}^{m},$

for some ${\omega}\in{}^{m}$ and, that ${\omega}$ is rationally
independent212121i.e. for any $k\in{{\mathbb{Z}}^{d}}\setminus\\{0\\}$,
${\omega}\cdot k\not=0$. if $m\geq 2$ and
${\omega}\in\real\setminus{\mathbb{Q}}$ if $m=1$. Then, the torus
$\phi({\mathbb{T}}^{m})$ is isotropic i.e. $\mathit{i}^{*}\varpi\equiv 0$,
where $\mathit{i}\colon\phi({\mathbb{T}}^{m})\hookrightarrow\mathcal{M}$ is
the inclusion map.

##  C Computation of the “action–angle” coordinates for elliptic billiards

It is well–note that the caustic parameter ${\lambda}$, together with its
conjugate $t\coloneqq F\big{(}\phi,k_{j}(\sqrt{{\lambda}})\big{)}$ is an
action–angle coordinates for the billiard in the ellipse
${\Omega}_{j}$.222222See e.g. [CF88] or [KS18, Proposition 16]. We are only
left to show that, ${\lambda}=g_{j}(\phi,{\theta})$ with $g_{j}$ defined as in
(19) is indeed the parameter of the caustic corresponding to the billiard
trajectory $\\{f_{j}^{n}(h_{j,1}(\phi,{\theta})):\;n\in{\mathbb{Z}}\\}$
starting at $h_{j,1}(\phi,{\theta})$. Let232323For the sake of simplicity, we
will drop the subscripts $j$.
$q_{\lambda}(t)=(q_{\lambda}^{x}(t),q_{\lambda}^{y}(t))=(a{\,\rm
cn}(t;k_{\lambda}),b{\,\rm sn}(t;k_{\lambda}))=(a\cos\phi,b\sin\phi)$, where
${\,\rm cn}$ and ${\,\rm sn}$ are the Jacobi elliptic functions:

${\,\rm cn}(t;k_{\lambda})\coloneqq\cos({\,\rm
am}(t;k_{\lambda}))\qquad\mbox{and}\qquad{\,\rm
sn}(t;k_{\lambda})\coloneqq\sin({\,\rm am}(t;k_{\lambda})),$

where ${\,\rm am}(t;k)$, called the amplitude of $t$, is given by:
$t=F\big{(}\phi,k)\big{)}$ iff $\phi={\,\rm am}(t;k)$. Then,242424see [CF88]
or [KS18, Proposition 16]. the trajectory starting at $q_{\lambda}(t)$ and
tangent to the caustic

${\Omega}_{{\lambda}}\coloneqq\left\\{(x,y)\in{{}^{2}}\ :\
\frac{x^{2}}{a^{2}-{\lambda}^{2}}+\frac{y^{2}}{b^{2}-{\lambda}^{2}}=1\right\\}\;,\qquad{\lambda}\in[0,b),$

hits the boundary at $q_{\lambda}(t+{\delta}_{\lambda})$, where
${\delta}_{l}\eqqcolon{\delta}({\lambda})$ is given in Lemma 21. Then, the
point of tangency $p_{\lambda}(t)$ of the segement
$[q_{\lambda}(t),q_{\lambda}(t+{\delta}_{\lambda})]$ to the caustic
${\Omega}_{\lambda}$ is given by
$p_{\lambda}(t)=(\sqrt{a^{2}-{\lambda}^{2}}\cos{\varphi},\sqrt{b^{2}-{\lambda}^{2}}\sin{\varphi})$,
for some ${\varphi}\in[0,2\pi)$. Thus, the slope of the line
$(q_{\lambda}(t),q_{\lambda}(t+{\delta}_{\lambda}))$ is given by:

  * $\bullet$

$-\sqrt{\frac{b^{2}-{\lambda}^{2}}{a^{2}-{\lambda}^{2}}}\cot{\varphi}$, using
its tangency property to ${\Omega}_{\lambda}$ at $p_{\lambda}(t)$;

  * $\bullet$

$\frac{b\sin\phi-\sqrt{b^{2}-{\lambda}^{2}}\sin{\varphi}}{a\cos\phi-\sqrt{a^{2}-{\lambda}^{2}}\cos{\varphi}}$,
using its points $q_{\lambda}(t)$ and $p_{\lambda}(t)$;

  * $\bullet$

$\tan\left({\theta}-\arctan\frac{b\cos\phi}{a\sin\phi}\right)=\frac{a\tan\phi\tan{\theta}-b}{a\tan\phi+b\tan{\theta}}$,
using the fact that the slope of the tangent to ${\Omega}$ at $q_{\lambda}(t)$
is $-\frac{b\cos\phi}{a\sin\phi}$ and that the angle between this tangent and
the line $(q_{\lambda}(t),q_{\lambda}(t+{\delta}_{\lambda}))$ is ${\theta}$.

Hence,

$\displaystyle-\sqrt{\frac{b^{2}-{\lambda}^{2}}{a^{2}-{\lambda}^{2}}}\cot{\varphi}$
$\displaystyle=\frac{b\sin\phi-\sqrt{b^{2}-{\lambda}^{2}}\sin{\varphi}}{a\cos\phi-\sqrt{a^{2}-{\lambda}^{2}}\cos{\varphi}}\,,$
(D.1)
$\displaystyle-\sqrt{\frac{b^{2}-{\lambda}^{2}}{a^{2}-{\lambda}^{2}}}\cot{\varphi}$
$\displaystyle=\frac{a\tan\phi\tan{\theta}-b}{a\tan\phi+b\tan{\theta}}\,.$
(D.2)

Now, (D.1) is equivalent to:

$a\sqrt{b^{2}-{\lambda}^{2}}\cos\phi\cos{\varphi}+b\sqrt{a^{2}-{\lambda}^{2}}\sin\phi\sin{\varphi}=\sqrt{(a^{2}-{\lambda}^{2})(b^{2}-{\lambda}^{2})}\,.$
(D.3)

The relation (D.2) implies

$\displaystyle\frac{\cos^{2}{\varphi}}{(a^{2}-{\lambda}^{2})(a\tan\phi\tan{\theta}-b)^{2}}=\frac{\sin^{2}{\varphi}}{(b^{2}-{\lambda}^{2})(a\tan\phi+b\tan{\theta})^{2}}=$
$\displaystyle\qquad\qquad\qquad=\frac{\cos^{2}{\varphi}+\sin^{2}{\varphi}}{(a^{2}-{\lambda}^{2})(a\tan\phi\tan{\theta}-b)^{2}+(b^{2}-{\lambda}^{2})(a\tan\phi+b\tan{\theta})^{2}}$
$\displaystyle\qquad\qquad\qquad=\frac{1}{(a^{2}-{\lambda}^{2})(a\tan\phi\tan{\theta}-b)^{2}+(b^{2}-{\lambda}^{2})(a\tan\phi+b\tan{\theta})^{2}}$
(D.4)

so that

$\displaystyle\frac{\cos{\varphi}}{\sqrt{a^{2}-{\lambda}^{2}}(a\tan\phi\tan{\theta}-b)}\stackrel{{\scriptstyle{\rm(\ref{slop2})}}}{{=}}\frac{-\sin{\varphi}}{\sqrt{b^{2}-{\lambda}^{2}}(a\tan\phi+b\tan{\theta})}=$
$\displaystyle\qquad\quad\stackrel{{\scriptstyle{\rm(\ref{sloop1})}}}{{=}}\frac{\pm
1}{\sqrt{(a^{2}-{\lambda}^{2})(a\tan\phi\tan{\theta}-b)^{2}+(b^{2}-{\lambda}^{2})(a\tan\phi+b\tan{\theta})^{2}}}\,.$
(D.5)

Thus, plugging into (D.3) the expressions of $\cos{\varphi}$ and
$\sin{\varphi}$ obtained from (D.5) and then squaring both sides of the
relation obtained yields ${\lambda}=g(\phi,{\theta})$, which completes the
computation.

## References

  * [BHS09] Hendrik W Broer, George B Huitema, and Mikhail B Sevryuk.  Quasi-periodic motions in families of dynamical systems: order amidst chaos.  Springer, 2009.
  * [Bia93] Misha Bialy.  Convex billiards and a theorem by e. hopf.  Mathematische Zeitschrift, 214(1):147–154, 1993.
  * [Bir13] George D Birkhoff.  Proof of poincaré’s geometric theorem.  Transactions of the American Mathematical Society, pages 14–22, 1913.
  * [Bir20] George D Birkhoff.  On the periodic motions of dynamical systems.  In Hamiltonian Dynamical Systems, pages 154–174. CRC Press, 2020\.
  * [CF88] Shau-Jin Chang and Richard Friedberg.  Elliptical billiards and Poncelet’s theorem.  Journal of mathematical physics, 29(7):1537–1550, 1988.
  * [Hal77] Benjamin Halpern.  Strange billiard tables.  Transactions of the American mathematical society, 232:297–305, 1977.
  * [Her89] Michael R Herman.  Inégalités a priori pour des tores lagrangiens invariants par des difféomorphismes symplectiques.  Publications Mathématiques de l’IHÉS, 70:47–101, 1989.
  * [Inc] Wolfram Research, Inc.  Mathematica, Version 11.3.
  * [Kal03] Vadim Kaloshin.  Geometric proofs of Mather’s connecting and accelerating theorems.  Topics in dynamics and ergodic theory, 310:81–106, 2003.
  * [KP90] Valery Kovachev and Georgi Popov.  Invariant tori for the billiard ball map.  Transactions of the American Mathematical Society, 317(1):45–81, 1990.
  * [KS18] Vadim Kaloshin and Alfonso Sorrentino.  On the local Birkhoff’s conjecture for convex billiards.  Annals of Mathematics, 188(1):315–380, 2018.
  * [Mat90] John N Mather.  Differentiability of the minimal average action as a function of the rotation number.  Boletim da Sociedade Brasileira de Matematica-Bulletin/Brazilian Mathematical Society, 21(1):59–70, 1990.
  * [MF94] John N Mather and Giovanni Forni.  Action minimizing orbits in hamiltomian systems.  In Transition to chaos in classical and quantum mechanics, pages 92–186. Springer, 1994.
  * [MM82] Shahla Marvizi and Richard Melrose.  Spectral invariants of convex planar regions.  Journal of differential geometry, 17(3):475–502, 1982.
  * [Sib04] Karl Friedrich Siburg.  The principle of least action in geometry and dynamics.  Number 1844. Springer Science & Business Media, 2004.
  * [Sie97] Martin Sieber.  Semiclassical transition from an elliptical to an oval billiard.  Journal of Physics A: Mathematical and General, 30(13):4563, 1997\.
  * [Sor15] Alfonso Sorrentino.  Computing Mather’s $\beta$–function for Birkhoff billiards.  Discrete & Continuous Dynamical Systems-A, 35(10):5055, 2015.
  * [Tab05] Serge Tabachnikov.  Geometry and billiards, volume 30.  American Mathematical Soc., 2005.