Analytical methods in Celestial Mechanics: satellites' stability and galactic billiards

In this paper, two models of interest for Celestial Mechanics are presented and analysed, using both analytic and numerical techniques, from the point of view of the possible presence of regular and/or chaotic motion, as well as the stability of the considered orbits. The first model, presented in a Hamiltonian formalism, can be used to describe the motion of a satellite around the Earth, taking into account both the non-spherical shape of our planet and the third-body gravitational influence of Sun and Moon. Using semi-analytical techniques coming from Normal Form and Nekhoroshev theories it is possible to provide stability estimates for the orbital elements of its geocentric motion. The second dynamical system presented can be used as a simplified model to describe the motion of a particle in an elliptic galaxy having a central massive core, and is constructed as a refraction billiard where an inner dynamics, induced by a Keplerian potential, is coupled with an external one, where a harmonic oscillator-type potential is considered. The investigation of the dynamics is carried on by using tools of ODEs' theory and is focused on studying the trajectories' properties in terms of periodicity, stability and, possibly, chaoticity.


Introduction
Phenomena involving the motion of Celestial bodies, either on a planetary or a galactic scale, are often characterised by a complex behaviour, whose accurate study requires the use of different tools, like numerical integration, analytical study, direct observations and much more ([4, 16, 9, 48, 2]).Analytical techniques represent, whenever applicable, useful strategies to study some of the main properties of the orbits in a gravitational system, especially in terms of long-term dynamics, providing results which are rigorous, as they follow form precise mathematical statements, and often general, in the sense that they potentially hold for a large set of trajectories (equivalently, of initial conditions).Along with purely analytical techniques, in some circumstances a mixed approach which includes numerics as well is possible: this is what happens, for example, when theoretical results are compared with simulations and observations, or in the case of semi-analytical approaches.Generally speaking, such expression refers to a class of methods where rigorous mathematical theorems are applied to numerically computed quantities (for example, within a Hamiltonian framework, to functions expressed through a truncated Taylor expansion, cfr.Section 2) This paper aims to illustrate the potential of such techniques, either analytic or semi-analytic, presenting the dynamical investigation of two models describing the geocentric motion of an object around the Earth and the trajectories of a body inside an elliptic galaxy with a massive core.In both cases, the motion of our test particle is influenced by the gravitational attraction of a variety of different mass distributions, depending on the model itself: as expected, the resulting dynamics is quite complex, and our main objective is to study its properties for long (possibly infinite) time scales.The issue of the long-term stability in geocentric motions is the core topic of Section 2, where a point-mass particle subjected to the attraction of the (non-spherical) Earth, Sun and Moon is taken into account.In general, the main question we try to answer is for how long it is possible to control the variation in the orbital elements (semimajor axis a, eccentricity e, inclination i) of our object, considering different initial conditions and, in particular, for different altitude regimes (we will use the classical distinction between NEO, MEO and GEO distances).Producing stability estimates for orbiting bodies at different distances from our planet's surface is a key problem in Celestial Mechanics, which finds its application to many different cases of practical interest.In particular, this problem is crucial when dealing with the wide and varied class of the object orbiting around the Earth, from satellites to microscopical space debris: in view of their large overall number and the collision hazard (cfr.[1]), the effort in predicting as well as possible their long-time behaviour has involved a remarkable community of mathematicians and astronomers (see for example [20,21,23,58,64], and, for a survey on the possible methods, [13]).Following the vast literature on the subject, the satellites' dynamics is formalised within a Hamiltonian setting via the so-called lunisolar Hamiltonian where r, ṙ are the position and velocity vector of our satellite in a suitable reference system, while the three potential parts refers to Earth's geopotential up to J 2 -term (see Section 2.1 and [46] for details) and Sun's and Moon's gravitational attractions treated as third-body perturbations.
As already anticipated, the techniques we used to produce stability estimates fall into the category of the semi-analytical methods, since H, opportunely treated, and expressed in its secular form (namely, averaged over the fast motions, see Section 2.1) can be written as a truncated Taylor expansion whose coefficients are computed numerically.From a mathematical point of view, we propose two different methods to produce stability estimates, holding in different regimes and based on different analytical results.With the first strategy, we provide stability estimates for the quantity I = √ µ E a √ 1 − e 2 (1 − cos i) in quasi circular and quasi equatorial orbits, computing an upper bound for the time up to which the variation of such quantity remains bounded within a certain range.The technique used to produce such estimates (see also [65]) is based on the application of a normal form algorithm: in short, and postponing the complete description to Section 2.2, we use canonical transformations in action-angle coordinates to reduce the Hamiltonian describing the satellites' motion to the form (1.1) where h 0 admits I as a first integral and the size of h 1 is so small that the overall dynamics can be considered a perturbation of the one induced by h 0 .The stability of I along the trajectories can be then deduced from the size of h 1 .The stability results, holding for small values of the eccentricity and inclination and five different values for the semimajor axis, corresponding to NEO, MEO, GEO region and beyond, are shown in Section 2.2, and precisely in Table 1.In short, one can say that the numerically computed stability times are extremely long, of the order of 10 4 years even in the worst case, corresponding to the farthest objects, although a worsening, due to the strengthening of the influence of Sun's and Moon's attraction, is evident beyond GEO region.
As for the second method, which is the subject of the analysis carried on in Section 2.3, it is based on Nekhoroshev theorem on exponential stability estimates (see [57]), and allows to cover a larger domain in eccentricity and inclination for satellites in MEO, and in particular for distances (in terms of semimajor axis) between 11 000 km and 19 000 km. Nekhoroshev theorem has already been used in some problems coming from Celestial Mechanics, like for example in the model of the Trojan asteroids ( [38]) and in the three-body problem ( [17,19]), and applies again to the case of a quasi-integrable Hamiltonian.Given an Hamiltonian function as in Eq.
(1. 1), where all the actions are first integrals for the unperturbed dynamics, under suitable conditions it provides stability times in the complete model which are exponentially long in the perturbation's size; the hypotheses required involve suitable nondegeneracy conditions on the unperturbed Hamiltonian h 0 , as well as a smallness condition on the size of h 1 .
In the present work, a nonresonant version of the theorem, which does not apply close to the secular geolunisolar resonances (see [12]), has been used, although a more complete analysis, which covers a wider regime, is possible provided a rigorous analysis of the geometry of resonances of the geolunisolar problem is carried on.The results in terms of stability times, for different values of a, e and i, are presented in Section 2.3, Figure 2.1: they are particularly good for low altitudes, and tend to worsen for increasing values of a.This phenomenon, that partially depends on the algorithm used to produce our stability estimates, will find an heuristic explanation in Section 2.3.
The second model taken into consideration is a simplified model that can be used to carry on a preliminary analysis on the motion of a particle in an elliptic galaxy having a central mass (a Black Hole or, in general, a massive core).This kind of motion, especially under the influence of super-massive bodies such as Black Holes, is particularly complex, and having a rigorous and reliable model to describe it would require to take into consideration the anisotropies in the mass distribution inside the galaxy, as well as relativistic effects.Situations of galaxies presenting Black Holes at their centers are quite common in actual galactic systems (see the review [34]), and it is quite natural, for anyone working in Celestial Mechanics, to ask how the presence of such a large central mass affects the dynamics, as well as whether it could lead to chaotic phenomena.The study we propose Section 3 has been inspired by the work [29], where the central body acts as a Keplerian center and the elliptic distribution of mass produce a harmonic oscillator-type potential.The superimposition of such two potentials leads to the establishment of two different regimes: whenever our test particle is close to the central body, the Keplerian attracting force of the latter is much more intense than the one of the overall galaxy, while the contrary happens whenever the particle is sufficiently far from the Black Hole.When the galaxy's mass distribution is an uniform ellipsoid, one can model its gravitational attraction via a harmonic oscillator-type potential (see [22]), where the frequencies over the three axes of the oscillations depend on the three semiaxes, while, ignoring possible relativistic effects, the potential of the central mass is a classical Keplerian one.In [29], the investigation of the model is carried on by means of a mixed analytical and numerical approach, and evidences of chaotic behaviour based on estimates of the corresponding Lyapunov exponents (see [43]) are shown.In our work (see also [27,28,7]), we propose a rigorous analysis which, although substantiated by numerical evidences, relies on a purely analytical approach: the price to pay is the necessity to introduce a simpler model, where the superposition of the two potential does not occur anymore.As in [29], we suppose that the distribution of matter in the ellipsoid is constant (except for the central body), and we build the model in such a way that our test particle is either attracted by the central mass or by the overall elliptic mass distribution.In practice, we divide the space into two regions, in each of which one of the two limit regimes identified by [29] occurs.In the inner region, representing the region of influence of the central mass, only its (Keplerian) attraction is considered; on the contrary, in a second region, exterior, the particle moves only Examples of orbits of refraction galactic billiards.The orbit goes inside and outside the domain, being deflected at every passage through the interface.Left: three-periodic trajectory.Right: quasi-periodic trajectory (figure taken from [28]).
under the influence of an isotropic harmonic oscillator.On the interface that separates these two regions the potential governing the particle's motion is generally discontinuous: to treat such discontinuity, we suppose that every time the particle hits the interface it undergoes a refraction, which deflects its velocity by a quantity that depends on the potentials' values in the transition point, as well as the hitting angle.Such refraction law, which in practice is a generalisation of classical Snell's law for light rays, can be interpreted as a limit case for a smooth passage from one potential to the other, where the intermediate region in which the two potentials are superimposed shrinks more and more (from a practical point of view, the paper [29] also provides estimates from the distance from the central mass at which the two potentials produce comparable forces).At this stage, we restrict our analysis on a planar system, considering one of the three invariant planes of the system identified by the ellipsoid's axes; this assumption will be removed in the next future ( [24]), where the three dimensional case will be considered.From a mathematical point of view, we construct our model by relying on the well established theory of mathematical billiards (see for example [66] for an extensive survey on the classical theory), and constructing the so-called galactic refraction billiard.In the classical case of mathematical billiards, a free particle moves inside a regular domain, following straight lines and bouncing against the boundary with an elastic reflection; deriving the properties of the particle's motion (equilibrium trajectories, existence of periodic orbits, chaotic regimes etc.) is a highly nontrivial problem, which has involved a wide community of mathematicians for at least one century.Within this framework, our billiard can be considered as a variation of the classical case, where the inner mass' domain of influence represents our billiard table, although two important differences have to be highlighted.First of all, the particle can exit from the domain, interacting with its boundary not with a simple reflection, but rather with a refraction that deflects its velocity vector.On the other hand, the presence of the inner and outer gravitational interactions leads to the appearance of two (outer and inner) non constant potential, so that the particle moves through curved geodesics instead of straight lines.Other examples of billiards with potentials, both in the reflective case and in the case of a coupled dynamics, are given in [36,50]; a remarkable example is given by Kepler billiards (see [59,67] and references therein), which, as we will see in Section 3.4, present strong analogies with our model.The results summarised in the present paper regard different aspects of the dynamics of the refraction galactic billiard for different domain's shapes and energy regimes, and have been achieved by using a wide class of tools coming from nonlinear analysis and the general theory of dynamical systems, as well as, sometimes, substantiated by numerical simulations.They will be presented into three main subgroups; first of all, as natural while dealing with a new dynamical system, the problem of existence and stability of equilibrium trajectories is considered.This is the topic of Section 3.2, where a particular class of equilibrium orbits, called homothetic and composed by straight lines, is considered.Such trajectories always exist when our domain is convex and smooth, and their linear stability can be studied by relying on the formalism of classical billiards and variational methods.In this framework, nontrivial bifurcation phenomena, occurring for non-circular domains with sufficiently smooth boundary, are shown, both from an analytical and a numerical point of view.After this preliminary (and local) analysis on the trajectory orbits, the investigation becomes more global in Section 3.3: here, the problem of the existence of periodic and quasi-periodic trajectories (see Figure 1.1) is treated.We work in a quasi-integrable regime, considering domains whose boundary is close to a circumference, and, as a consequence, can be treated with the powerful tools of perturbation theory (see also [13], where such concepts are explained in a slightly different framework).In particular, we shall make use of KAM theorem (see [55]), Poincaré-Birkhoff and Aubry-Mather theories (see [39]) to prove that, whenever our domain's boundary is smooth enough and sufficiently close to a circle, then there exists orbits with any rotation number within a certain range (see Theorem 3.7 and Eq.(3.12) for the formal definition of rotation number).We stress that our case is not the first application of KAM, Aubry-Mather and Poincaré-Birkhoff theories in problems coming from Celestial Mechanics: examples are [14,10].The landing point of the analysis of the galactic billiards' dynamics, at least in the regime here presented, is included in Section 3.4, where we address the problem of the possible chaoticity of the system.In our specific case, evidences of chaotic behaviour are included in both [29] and the numerical simulations presented in [27] and reproduced in Figure 3.9: these two elements motivated the prosecution of the study in this direction, trying to formally prove the chaoticity of the model.The final result resumes in the detection of a simple geometric condition on the domain's shape, called admissibility, that ensures the existence of a topologically chaotic subsystem of the galactic refraction billiard for large enough inner energies (see Theorem 3.9).Roughly speaking, and postponing the rigorous definition to Section 3.4 (and in particular Definition 3.8), we say that a domain with smooth boundary is admissible whenever there exist two segments from the Keplerian mass which are orthogonal to the boundary, not antipodal with respect to the origin and such that they are nondegenerate (that is, the hitting point is a strict maximum or minimum of the function distance from the mass restricted to the domain's boundary), see Section 3.4, Figure 3.10.In practice, admissibility acts as a sufficient condition that, through a particular construction called symbolic dynamics, ensures that, up to restricting to a subset of all the possible initial conditions and choosing a sufficiently large inner energy, our galactic billiard is chaotic.Also in this case, we stress that our work can be considered as a part of a vast literature, whose aim is to investigate and detect, with different techniques, chaotic systems in Celestial Mechanics, both with a rigorous analytical approach (see for example [8,2,3]), and with a more numerical point of view [40,35].

Hamiltonian methods for satellites' stability estimates
The current Section summarises the results of [26,15] regarding the long-term stability for bodies orbiting around the Earth, considering, in a Hamiltonian setting, the gravitational attraction of our planet, Sun and Moon.Section 2.1 describes the Hamiltonian model taken into consideration, including the set of action-angle variables used, with particular attention on their physical meaning in terms of orbital elements.Section 2.2 resumes the main ideas behind normal form theory, proposing then the application of such approach to our model to produce stability estimates for eccentricity and inclination, locked in the quasi-integral I = √ µ E a √ 1 − e2 (1 − cos i), for quasi circular and quasi equatorial orbits.Section 2.3 widen the set of the considered initial conditions to more inclined and eccentric orbits within MEO distances: in this case, an approach based on the application of Nekhoroshev theorem is taken into account.Finally, Section 2.4 presents some final considerations on the results obtained, comparing the two approaches both in terms of the numerical outcome and theoretical consequences.
2.1.The Hamiltonian model.To construct the Hamiltonian function related to the geolunisolar model, let us start by considering a geocentric reference frame in the space, with coordinate axes x, y and z, where the (x, y)−plane corresponds to the Earth's equatorial one and the x−axis points towards the line of the equinox.In such framework, the geolunisolar Hamiltonian referred to a point-mass particle of coordinate vector r = (x, y, z) can be expressed by1 where the potential terms are given as follows: • the term H E is the Earth's gravitational potential which takes into account the nonspherical shape of our planet; it can be expressed as an expansion in spherical harmonics, as described in [45].In the current model, such expansion is truncated up to the J 2term 2 , giving rise to an expression of the form where R E = 6378.14km and µ E = GM E = 1.52984 × 10 9 R 3 E /yr 2 are respectively the Earth's radius and mass parameter and J 2 = −1082.6261× 10 −6 is a dimensionless parameter; • the terms H S and H M refer to the gravitational attraction of Sun and Moon, whose motion in the geocentric reference frame is given respectively by the time-dependent position vectors r S (t) = (x S (t), y S (t), z S (t)) and r M (t) = (x M (t), y M (t), z M (t)).More precisely, one has where again µ S and µ M are respectively the mass parameters of Sun and Moon.As for the analytic expression of r S and r M , one has that both bodies moving around the Earth describing ellipses: the orbital parameters (inclination i 0 , semimajor axis a and eccentricity e) of Sun are i 0S = 23.43 • , a S = 1.469 × 10 8 km and e S = 0.0167, while for the Moon one has i 0M = i 0S , a M = 38 4748 km and e M = 0.065.From a dynamical point of view, assuming that the Moon lies on the ecliptic plane corresponds to neglect the precession of the Lunar node: as will be observed later in Section 2.3, this assumption will have important effects close to the so-called secular lunisolar resonances (see also [12]).
Since the motion of our point-mass particle is a geocentric trajectory, it is convenient to express the Hamiltonian 2.1 in terms of the particle's orbital elements; such change of variables is performed by expressing, as in [25] and [56], the coordinates x, y, z (resp.the components x S\M , y S\M , z S\M of r S\M ) in terms of orbital elements (a, e, i, M, ω, Ω) (resp.a S\M , e S\M , i S\M , M S\M , ω S\M and Ω S\M ), where a, e and i denote respectively the orbit's semimajor axis, eccentricity and inclination, while the angles M, ω and Ω are the mean anomaly, the argument of the perigee and the longitude of the nodes.The resulting Hamiltonian, which will still be called H, is a function of (a, e, i, M, ω, Ω), where where the time dependence is expressed by the mean anomalies M S and M M of Sun and Moon.When one in interested to the satellite's long term dynamics, the Hamiltonian H can be further simplified, removing the dependence on the time, by considering an averaging process over the fast angles of the problem (namely, the three mean anomalies): the result of this averaging is the secular geolunisolar Hamiltonian (2.2) The Hamiltonian H sec is the starting point to obtain long-term stability estimates for the secular geolunisolar model, either with normalization techniques, as in Section 2.2, or through the application of stronger results, such as Nekhoroshev Theorem, as in the case of Section 2.3.
In order to carry on such investigation, one needs to express the above Hamiltonian in terms of action-angle variables, such as the so-called modified Delaunay ones (see [54]), whose relation with the orbital elements is given by Note that, in terms of these new variables, the averaging performed above corresponds to the elimination of the fast angle λ, and, subsequently, to take the first action L (and then the semimajor axis) as a constant, which we call L * = √ µ E a * .As for the eccentricity and inclination, the presence of Sun and Moon on the ecliptic forces the existence of a circular, non-equatorial equilibrium orbit, whose inclination i (eq) depends on a * through the relation The equilibrium points (e (eq) , i (eq) ), with e (eq) = 0, are traditionally called the forced elements of the secular model, while the plane with inclination i (eq) is the Laplace plane; more rigorous estimates on the value of i (eq) and its behaviour as a function of the distance can be found in [62].
The stability estimates produced in this work are obtained by means of a semi-analytical approach, namely, the application on rigorous analytical results on Hamiltonian computed numerically by means of the software Mathematica © : for this reason, in Sections 2.2 and 2.3 we shall make use of a truncated expression of H sec , whose truncation order will be specified case by case.

Stability estimates through normal forms.
The first technique we propose to estimate the stability of the orbital elements in the secular geolunisolar model relies on the application of a normal form algorithm, and is similar to the one used in [65].Before passing to the actual computation of the stability time in the satellites' case, a brief general introduction of the normal form theory is in order (a more complete dissertation on the subject can be found in [25,32]).Let us start by taking a Hamiltonian function expressed in action-angle variables H(J, θ), where (J, θ) ∈ U × T n , n being the degrees of freedom of the system and U ⊂ R n open.The principal aim of a normalization algorithm is to find a close-to-identity canonical transformation Φ : (J, θ) → (J ′ , θ ′ ) such that the new Hamiltonian H ′ = H • Φ −1 takes the form where: • Z(J ′ , θ ′ ) is the so-called normal part, and has some desired property as, for example, the presence of first integrals of the motion; • R(J ′ , θ ′ ) is the remainder : in a suitable functional norm ∥•∥, it is such that ∥R∥ ≪ ∥Z∥.If the remainder's size is sufficiently small with respect to the normal part's ones, the overall dynamics under H ′ (and then under H) can be considered as a small perturbation of the one induced by Z.As an example (which will be precisely our case), if Z admits some integrals of the motion, such quantities are quasi-constant for the whole H ′ .The transformation Φ can be found by means of the Lie series technique: its construction algorithm, which depends on the properties of the normal part we seek, is here omitted, and can be found in [32].In this Section, a normalization algorithm is used to produce stability estimates for the eccentricity and inclination for orbits close to the equilibrium one, which has orbital parameters (a * , e (eq) , i (eq) ) (see Section 2.1).As a preliminary step for this analysis, it is convenient to consider a set of Delaunay coordinates which are centered around the equilibrium, performing the change of coordinates (2.4) where P (eq) = √ µ E a * 1 − 1 − (e (eq) ) 2 = 0 and Q (eq) √ µ E a * 1 − (e (eq) ) 2 1 − cos i (eq) = √ µ E a * 1 − cos i (eq) .By means of a Taylor expansion, the Hamiltonian can be then written as a trigonometric polynomial in the square roots of the actions as (2.5) By construction and from Eq.( 2.3), one has that, for quasi-circular orbits close to the Laplace plane (2.6) where (e, i) have to be intended as the differences with respect to the forced values (0, i (eq) ); this implies that, in the expansion (2.5), the s−th term in the sum is of total order s in eccentricity and inclination.For computational reasons, in the following estimates the series in Eq.(2.5) is truncated up to order N = 15.The rigorous procedure to obtain H(I 1 , I 2 , ϕ 1 , ϕ 2 ) is described in [25], where one can also observe that the first order frequencies ν 1 and ν 2 are nearly equal : this fact, which implies a 1 : 1 resonance between the conjugate angles ϕ 1 and ϕ 2 , will be crucial in the normalization procedure.
Once the Hamiltonian is in the form of Eq. ( 2.5), one can proceed with the normalization: in this case, it consists in finding a change of coordinates which makes the normal part depending only on the resonant angle ϕ 1 − ϕ 2 , namely, a near-identity canonical transformation Φ such that the new Hamiltonian (which, with an abuse of notation, will be still called is given by the sum (2.7) From a practical point of view, this result is achieved by applying to the initial Hamiltonian a sequence of transformations in the form of Lie series, aiming to remove the dependence on the angles, except for the resonant combination ϕ 1 − ϕ 2 , from all the terms of the series (2.5) up to order M = 12.The choice of the normalization order M is of particular importance to obtain optimal estimates: for more details on that, see [33].
Hamiltonians in the form of (2.7) are usually said to be in resonant normal form: here, the normal part is composed by the secular term Z sec , which does not depend on the angles, and the resonant one, that depends on the actions as well as on the resonant combination ϕ 1 − ϕ 2 ; as for the remainder, it is by construction of order M in the square roots of the actions (namely, recalling Eq.(2.6), in eccentricity and inclination3 ).
It is easy to prove that the quantity 2 is an integral of the motion for the dynamics induced by the sole normal part Z sec + Z res : the conservation of such quantity, which is equal to the vertical component of the satellite's angular momentum, determines a locking between eccentricity and inclination, which can undergo only changes which keep constant the value of I 1 +I 2 .This fact, also known as Lidov-Kozai effect (see [49,52]), is common in many model of Celestial Mechanics which present resonance phenomena.For the overall dynamics induced by (2.7), the quantity I 1 + I 2 it is not constant anymore; nevertheless, if the remainder's norm is sufficiently small, is can be considered as quasi-constant, and it is possible to obtain stability estimates (namely, an upper bound for the time up to which it is bounded in a certain neighborhood around the initial values of e and i) by measuring the size of R in a suitable functional norm.More precisely, let us fix a domain D ⊂ R 2 around the forced values for eccentricity and inclination (0, i (eq) ), and, given a function f (e, i, Our final objective is to evaluate the variation of I 1 + I 2 along the trajectories induced by the normalized Hamiltonian in (2.7): to this aim, let us recall the relation where the notation {•, •} denotes the Poisson brackets (see [37]).Being I 1 + I 2 a first integral for the normal part, {I 1 + I 2 , Z sec + Z res } = 0, and then, given any (ê, î, φ1 , φ2 ) ∈ D × T 2 , one has (2.9) Let us now suppose that at the time t = 0 the quantity I 1 + I 2 has value I 0 1 + I 0 2 , corresponding to eccentricity and inclination (e 0 , i 0 ) ∈ D, and consider its time evolution over t.Suppose now to fix Γ > 0 as the maximal variation allowed over a certain time for (I 1 + I 2 )(t): applying the mean value theorem, it is possible to compute an upper bound for the time T such that, for any t ≤ T , More precisely, from Eq.(2.9) one has that .
The upper bound T = Γ/∥{I 1 + I 2 , R}∥ D,∞ is the stability time we seek: it depends of course on the maximal variation allowed Γ, as well as on the amplitude of the domain D in eccentricity and inclination we want to analyse.It is clear that the value of T increases with Γ and by taking smaller domains around the forced elements; moreover, it has a dependence on the reference value of the semimajor axis a * .
For computational reasons, to produce the numerical estimates on T the domain D is set to be D = {(e, i) ∈ [0, 0.1] × [0 rad, 0.1 rad]}, while Γ will depend on a * through the relation additionally, the sup norm in (2.8) is replaced with an alternative functional norm based on majorization (see the details in [25]).
It is clear that the whole stability estimate process depends crucially on the semimajor axis, which in the secular geolunisolar model is a constant parameter; for this reason, in the numerical estimates we distinguished within five different cases, which cover many different regimes (for the sake of clarity, they will be given in terms of the sum of the altitude with the Earth's radius): • a (1) * = 3 000 km + R E , corresponding to an orbit just above the atmosphere; • a (2)   * = 20 000 km + R E , located in the MEO region; • a (3)   * = 35 786 km + R E , which corresponds to the altitude of GEO orbits; • a (4)   * = 50 000 km + R E , corresponding to far object; • a (5)   * = 100 000 km + R E , which is the distance of objects very far from the Earth, where the influence of Sun and Moon is particularly strong.Table 1 shows the stability times obtained for these values of the semimajor axis; as one can easily notice, though particularly long, the time T decreases with the altitude, with a significant worsening beyond GEO distance.These results, obtained numerically, are consistent with the theory: for small values of the semimajor axis the secular geolunisolar model can be well approximated by the secular J 2 model (namely, the model in which only the geopotential up to the term J 2 is considered,

Semimajor axis Stability time T (years)
a , obtained via normalisation method, for quasi circular and quasi equatorial orbits of the geolunisolar model and five different reference values of the semimajor axis.The data are taken from [26] .averaged over the mean anomaly), which is integrable; on the other hand, going farther from the Earth's surface, the influence of Sun and Moon gets stronger and stronger, leading to a perturbation that produces instability in the model.

Exponential stability estimates through Nekhoroshev Theorem.
In addition to be used as in Section 2.2 to produce stability estimates based on the mean value theorem, normal form algorithms are an essential preliminary tool (as well as a proving strategy) to apply the celebrated Nekhoroshev theorem (see [57,60]), here presented in its nonresonant form.In general (we will be more precise in Theorem 2.1 for the specific case of the nonresonant regime), such theorem can be applied to quasi-integrable Hamiltonians of the form where h 0 depends only on the actions and h 1 depends on the angles as well.As a consequence of Hamilton's equations, the dynamics induced by h 0 has the actions as first integrals of the motion, namely, J(t) = J 0 for any t ≥ 0, J 0 being their initial value.Under some suitable nondegeneracy condition on h 0 and provided that the perturbative function h 1 is small enough, it is possible to estimate the stability time of the actions under the dynamics induced by the whole H: in particular, it is possible to find an open set around J 0 where the actions are bounded for a time which is exponentially long in the inverse of the perturbation's norm.
In the most general formulation of Nekhoroshev theorem, the nondegeneracy hypothesis required on h 0 , called steepness condition, resumes essentially in asking for a quantitative transversality condition for the gradient ∇h 0 ; here, we will rely on a simpler nonresonance hypothesis, based on the non-commensurability of the coefficients of the actions at first order, which can be easily verified numerically.As we will see while presenting the numerical results (see Figure 2.1), the application of this simpler version of the theorem implies a cost in terms of the region of the (a, e, i)−space where our estimates hold: nevertheless, the stability times obtained are particularly good in a strip of the MEO region and in a nonresonant regime, being comparable with the satellites' average orbital lifetime; a finer analysis, considering the geometry of the resonances in the geolunisolar problem, is anyway possible.
To apply the nonresonant version of Nekhoroshev theorem to our geolunisolar case, it is necessary to put the Hamiltonian (2.2) in the form of a sum of an integrable term and a perturbation, as in Eq. (2.10): to this aim, we will rely again on a normal form algorithm.
Hamiltonian preparation.Let us start again from the secular geolunisolar Hamiltonian as presented in Eq.(2.2).While in Section 2.2 we focused our investigation in a small neighborhood (in eccentricity and inclination) of the forced elements (0, i (eq) ), here we aim to provide stability times holding for values of the orbital parameters which are not necessarily small.We will then produce a sequence of Hamiltonian functions, each of which is obtained by expanding H sec around a triplet of reference values (a * , e * , i * ) in a grid covering the set [11 000 km, 20 000 km]× [0, 0.5] × [0 • , 90 • ]; for each Hamiltonian of such sequence, we will follow a numerical procedure, described below, to provide stability estimates holding in a neighborhood of the corresponding reference values (e * , i * ) (remember that, in the secular geolunisolar problem, the semimajor axis a * is a priori constant for any forward time).
In practice, once fixed (a * , e * , i * ), one can perform a translation in the actions analogous to the one presented in Eq. (2.4) to obtain an expansion of the form (2.11) where the j−th term of the sum is given by Note that the expansion in Eq. (2.11) is the analogous of Eq.(2.5) in Section 2.2, although in this case the exponential form has been chosen.The explicit expressions of ω 1 and ω 2 , as well as of the coefficient of a (j) and b (j) at first and second order, can be found in [15].By computing ω 1 and ω 2 numerically, it is possible to observe that there are particular values of the reference inclination i * for which they are commensurable: these are the inclinations of the so-called secular resonances for the geolunisolar problem (see for example [12]), which will play a fundamental role in the upcoming stability analysis.As in Section 2.2, for computational reason the sum in Eq.(2.11) has been truncated, up to order N = 12.The next step towards producing stability estimates via Nekhoroshev theorem consists in normalising the Hamiltonian in Eq. (2.11), namely, using canonical transformations to obtain an expression as in Eq.(2.10),where h 0 contains only angle-independent terms and the size of h 1 can be controlled with a suitable norm.Let us start by considering the non-normalised sum in Eq. (2.11), and suppose to split it into the form where h0 contains only the angle-independent terms in (2.11) and h1 contains all the others.If we suppose that the action values are bounded, is clear that, from j = 4 on, the size of the angle-dependent summands decrease quadratically with the action's bound; on the other hand, the purely trigonometric terms, as well as angle-dependent ones which are linear in the actions, are harder to control 4 .The normalisation algorithm performed in this case aims precisely to the elimination of such terms up to a certain order M (in the actual computation, M is set equal to 6) via a sequence of suitable Lie series transformations.The complete algorithm, whose extended description can be found in [15], leads finally to a new Hamiltonian (2.12) where, with an abuse of notation, the new action-angle variables are still called I i and ϕ i , i = 1, 2. In Eq.(2.12), the normal part, composed by the linear part plus Z, is composed by angle-independent terms plus other terms which could be angle-dependent but at least quadratic in the actions.As for the remainder term R, it could contain terms which depend on (ϕ 1 , ϕ 2 ) and are constant or linear in the actions; nevertheless, provided the normalization algorithm converges (namely, the size of the coefficients of the remainder decreases in the process), such terms are small with respect to Z.The convergence of the normalization is a crucial issue of the overall procedure, which depends heavily on the non-commensurability of the initial frequencies ω 1 and ω 2 ; furthermore, such convergence influences also the final value of the frequencies, denoted by ω1 and ω2 .The variation in such quantities is negligible whenever the normalisation converges.
The normalized Hamiltonian H norm is the starting point to obtain exponential stability estimates via non-resonant Nekhoroshev theorem, which we now recall in the version by Pöschel (see [60]), after some useful definitions.
Let us start by considering, in general, a Hamiltonian of the form which is assumed to be real analytic in (J, θ) ∈ A × T n , A ⊂ R n .Suppose also that the above Hamiltonian can be extended analytically to the set (2.13) where r 0 and s 0 are two positive real constants.As a last assumption, let us suppose that the Hessian matrix associated to h 0 is bounded in A r 0 , namely, that there exists a constant M > 0 such that, denoted with ∥ • ∥ o the operator norm induced by the Euclidean one on R 2 , sup Given now an analytic function expressed as we define the Cauchy norm of g |g| A,r 0 ,s 0 = sup then for every orbit of initial conditions (J 0 , θ 0 ) ∈ D × T n one has Once one has precise numerically computed values for all the quantities involved, one can use Theorem 2.1 to produce stability estimates for the actions; more precisely, in a non-resonant regime defined through the notion of α − K−nonresonance, one can find an open set in R n in which the actions are bounded for a time which is exponentially long in K.We stress that such result is to be intended as local, in the sense that it holds for initial values for the actions in a subset D of A. Numerical evidences (see the paragraph "Numerical results" below) show how the cut-off value K satisfies a relation of the type K ∼ (c 1 |h 1 | A,r 0 ,s 0 ) −c 2 , c 1 and c 2 being two positive constants.Such behaviour is consistent with theoretical results (see for example [60]), and allows to conclude that the stability time is exponentially long with respect to the perturbation's norm to some power.
Numerical results.To produce stability times though Theorem 2.1 for the secular geolunisolar model, it is necessary to set up an algorithm that, given reference values of the orbital elements, computes the quantities involved in the Theorem and finally, if the hypothesis (2.14) is satisfied, provides T stab as in Eq.(2.15).In practice, our algorithm develops into the below steps.
(1) We start by fixing the constants a, b, r 0 , s 0 , whose value has been established by trials and errors, and could be possibly tuned to obtain optimal estimates.In particular, we impose r 0 = s 0 = 0.1, a = 9/8 and b = 1/8 (the choice of a and b's values is the same one can find in [60]).Moreover, we fix a reference value of the semimajor axis a * , which, by virtue of the averaging process, is constant along every orbit.(2) Fixed the values (e * , i * ), we compute numerically the expansion (2.11), to arrive, after the normalization, to the form (2.12).We stress that the final values of the frequencies ω1 , ω2 , as well as the actual size of the remainder term R, depend heavily on the noncommensurability of ω 1 and ω 2 , namely, on the reference values (e * , i * ).This means that different values of eccentricity and inclination could lead to completely different outcomes in terms of normalization.The Hamiltonian can be now splitted into an integrable part h 0 (I 1 , I 2 ) containing only the angle-independent terms plus a perturbation h 1 (I 1 , I 2 , ϕ 1 , ϕ 2 ) which contains all the other terms.(3) We can now compute the quantities involved in Theorem 2.1: first of all, we define the actions' set A as 1 , I * 2 being the actions corresponding to the reference values; one can then define D r 0 ,s 0 as in Eq. (2.13) and As for the nonresonance parameters α, K, we search for their optimal values, provided condition (2.14) is satisfied, as follows: for every i = 1, . . ., 50 we compute At this point, one can compute |h 1 | A,r 0 ,s 0 and check whether there exists i ∈ {1, . . ., 50} such that |h 1 | A,r 0 ,s 0 ≤ ϵ * i : if it happens, then one can take K as the maximal i such that the condition is verified, α = α K and compute the stability time as in Eq.(2.15).On the other hand, since the sequence there is no hope for the theorem to be applied for the specific values (a * , e * , i * ) and any i ∈ {1, . . ., 50}: in this case, we impose K = 0.The color scale refers to the computed stability times (in years), while the white region correspond to the values of (e, i) where Theorem 2.1 can not be applied with the present algorithm.The red lines are in correspondence of the inclinations of the secular geolunisolar resonances.Data taken from [15].0.5 • in inclination.The color scale indicates the value of the stability time (in years) obtained, while the white region of such values of (e * , i * ) for which condition (2.14), using the proposed algorithm, does not hold.The red lines are put in correspondence of the known inclinationdependent resonances for the secular geolunisolar problem (see [12]).From the numerical results, it is evident as the domain where the Theorem can be applied shrinks manifestly with a * , and, concurrently, the estimates on the stability times get worse.Moreover, an evident influence of the resonances comes out, since, even in the best case (i.e. for a * = 11 000 km), white regions around the corresponding inclinations appear.The role played by the resonances in the overall procedure enters at two different levels: during the normalization process and, later, in the very application of Theorem 2.1.As for the first normalization, one can check from the explicit expression of the coordinate changes used to remove the "unwanted" terms from the normal part (see [15] for all the details) that linear combinations of ω 1 and ω 2 appear at the denominator: whenever the frequencies are resonant, such denominators (the so called small divisors) approach zero, leading to an explosion in the remainder R and, subsequently, in the size of h 1 .We refer again to [15] for a detailed analysis of the convergence of the first normalization, including results on the change in the frequencies' value during the process.
On the other hand, the simple fact that we are using a nonresonant version of Nekhoroshev theorem makes clear how having a commensurability relation at low order between ω1 and ω2 correspond to a value of α (and, as a consequence, of the threshold ϵ * ) drastically low, making nearly impossible for the norm of h 1 to remain below ϵ * .The effect of the distance on the worsening of the results has a more complex reason, which can be explained, roughly, by the following heuristic argument: it can be shown that, after the normalization, the remainder R contains purely trigonometric terms whose size is comparable to the one of Ca 5 * tan i * M , where C is a suitable constant and M is the normalization order, here put equal to 6.As a consequence, the size of these terms, which can not be controlled by taking a smaller domain in the actions, grows swiftly with a * and whenever i * approaches 90 • .We conclude the exploration of the numerical analysis of the stability problem by providing an example which shows the behaviour of the computed value of the cut-off value K with respect to the perturbation's norm |h 1 | A,r 0 ,s 0 .Figure 2.2 shows the LogLog plot of the values of K and |h 1 | A,r 0 ,s 0 for a * = 13 000 km, e * = 0.2 and inclinations in a mesh of [0 which is consistent with the expected theoretical results, and allows to conclude, as anticipated before, that the final estimates can be actually considered as exponentially long in the inverse of the perturbing function's norm.

2.4.
Further considerations and conclusions.The techniques used in Sections 2.2 and 2.3 are examples of how semi-analytical manipulations in a Hamiltonian framework could be used to gain information in the long-term dynamics of a body orbiting around the Earth under the influence of the latter, Sun and Moon.Other examples of this kind can be found in the literature (see for example [25], where also the case of the J 2 −model is taken into consideration, or [47,18,61]).The first method, inspired by the work of Giorgilli and Steichen in [65], essentially provides stability times which, though very long (see Table 1), are linear with respect to the perturbation's norm, and hold in quasi-circular orbits lying close to the Laplace plane.As for the second method, it produces estimates which are exponentially long with respect to the inverse of the perturbation's size, showing all the potential of the Nekhoroshev theorem (see [57]), which could be used also in higher dimensions; on the other hand, at present the domain in which the results are truly substantial is not particularly large (see Figure 2.1).Nevertheless, we stress that different strategies to obtain an initial normal form may overcome the convergence problem, and, most of all, that the above procedure is based on a nonresonant result: a finer analysis of the geometry of the resonances in the secular geolunisolar problem would allow to use Nekhoroshev theorem in its complete version, obtaining estimates valid in a resonant regime as well.
As for the model we chose to use, we stress that we are considering the influence of the geopotential only up the J 2 −term.The overall analysis can be refined by taking also further terms in (2.1), like for example the ones corresponding to J 2 2 , J 3 and J 4 .A comparison between our results and the ones one can obtain by considering this more complete model is presented at the end of [15], showing that, in the practical context on the satellites' motion, the stability times obtained for the two models, though different, are so long with respect to the average operational lifetime that any change does not really affect the validity of the estimates.

Regular and chaotic motions in Galactic Billiards
The current Section resumes the results contained in [27,28,7] on the analysis of the refraction galactic billiard (see Section 1), a model aiming to provide a simplified description of the motion of a particle in an ellipsoidal galaxy having a central super-massive core.Section 3.1 is intended as a description, in this framework, of the considered dynamical system, complete with the motivations that led us to choose a the considered potentials, as well as a refractive interface.The results obtained are divided into three subgroups: Section 3.2 is focused on the existence and linear stability of equilibrium trajectories for the model, and provides numerical and analytical evidences of bifurcation phenomena regarding a particular class of orbits.Section 3.3 extends the analysis to periodic and quasi-periodic trajectories, within a perturbative regime constructed by taking into consideration quasi-circular billiards.In Section 3.4, the problem of the possible arising of chaotic behaviours is taken into account, arriving to the detection of simple geometric conditions on the billiard's boundary that ensure the presence of a chaotic subsystem at high energies.
where E, ω, h, µ are positive constants representing respectively the energy and frequency of the outer harmonic oscillator, the difference in energy between inner and outer trajectories and the central body's mass parameter.Starting from initial conditions on the interface ∂D, the trajectories at zero energy induced by the inner potential are Keplerian hyperbolae, while the outer ones are elliptic harmonic arcs: with a broken geodesics technique (see [63]) we can construct complete trajectories in our system by patching together outer and inner arcs.The connection rule is given by the refraction Snell's law described in Figure 3.1, left: denoting with α I and α E respectively the angles of the inner and outer arcs connected at a point z ∈ ∂D with respect to the normal direction to ∂D in z, the following relation must be satisfied Geometrically, Eq.(3.1) translates in the conservation of the tangent component of the velocity after the transition.The choice of this kind of connection rule is based on different arguments: first of all, from a physical point of view, it can be seen as a generalisation for non-constant potentials and non-straight interfaces of the classical Snell's law for light rays.On the other hand, it has a rigorous and robust variational interpretation, which will be crucial in the whole forthcoming analysis.To explain it (see [27] for further details), let us consider a concatenation of an outer and a inner arc that connects two point on the boundary p 0 and p 1 , passing trough a transition point p (see Figure 3.1, right).It is possible to associate to any of the two arcs, denoted for the moment by z E (t) and z I (t), the corresponding Jacobi lengths where z E (0) = p 1 , z E (T E ) = z I (0) = p and z I (T I ) = p 1 .Under suitable conditions, it can be proved that the outer (resp.inner) trajectory under the potential V E (resp.V I ) arc connecting two points on the boundary is unique: as a consequence, the functions L E and L I depend only on the endpoints.The inner and outer Jacobi lengths can be combined to obtain the total Jacobi length of our concatenation, and, making use of this quantity, it is possible to state Snell's law in a variational way as follows: we say that the concatenation from p 0 to p 1 through p satisfies Snell's law at the transition point if and only if p is a critical point for the total Jacobi length of the concatenation itself, that is, Of course, an analogous reasoning applies whenever the transition is from inside to outside.
As customary in billiards theory, to study the two dimensional dynamics of the trajectories of the complete system it is possible to restrict ourselves to a discrete map which keeps track of the behaviour of a concatenation whenever it hits the boundary: this is the so-called first return map, which, starting from generic initial conditions on the boundary (position and velocity vector), summarises the behaviour of the generated trajectory after every concatenation of an α ′ ṽ z I (t; p, ṽ) First return map: starting from initial conditions (p 0 , v 0 ), determined by the one-dimensional parameters (ξ 0 , α 0 ), the trajectory is follow through an outer arc, a refraction from outside to inside, an inner arc and a refraction from inside to outside to find the final conditions (p 1 , v 1 ), defined by (ξ 1 , α 1 ).outer and subsequent inner arc.To be more precise, let us start by parametrising ∂D with a smooth, closed and simple curve γ : I → R 2 , ξ → γ(ξ), where I ⊂ R is a suitable interval.For the sake of simplicity and without loss of generality, we can suppose that γ is the arc length parametrisation of ∂D, so that | γ(ξ)| = 1 for any ξ ∈ I. Let us now take initial conditions on the boundary for an outer arc, (p 0 , v 0 ) ∈ ∂D × R 2 , such that v 0 points outside D and the energy conservation law for the outer problem is satisfied, that is, |v 0 | 2 /2 − E + ω 2 |p 0 | 2 /2 = 0 (see Figure 3.2).Such initial conditions are uniquely determined by a pair of one dimensional parameters (ξ 0 , α 0 ), with γ(ξ 0 ) = p 0 and α 0 the angle between v 0 and the outward-pointing normal unit vector to γ in ξ 0 .Once the initial conditions are fixed, we can consider the outer arc z E (•; p 0 , v 0 ) which is a solution of the outer Cauchy problem Since ∂D is bounded and z E is an elliptic arc, there exists a first return time for the outer dynamics; in other words, there exists T E > 0 such that z E (T E ; p 0 , v 0 ) ∈ ∂D and z E (t; p 0 , v 0 ) / ∈ D for any t ∈ (0, T E ).We can then consider ( ξ, α′ ) as the parameters that describe, using the same rationale as before, (z E (T E ; p 0 , v 0 ), z ′ E (T E ; p 0 , v 0 )).At this point, we have a trajectory that hits the boundary and must be refracted: following Eq.(3.1) with z = γ( ξ) and α E = α′ , one can find α I and the corresponding initial conditions for the outer arc (p, ṽ) (without giving the analytical formulae, we refer again to Figure 3.2), to define the inner arc z I (t; p, ṽ) as the solution of the ODE Note that Eq.(3.1) implies automatically that the initial velocity satisfies the inner energy equation, that is, |ṽ| 2 /2 − (E + h) − µ/|p| = 0. Again, given that z I is an unbounded Keplerian hyperbola, there exists a first return instant T I on ∂D: we can then take the outer arc's final conditions (z I (T I ; p, ṽ), z ′ I (T I ; p, ṽ)) and refract the inner velocity to obtain conditions (p 1 , v 1 ) such that it holds p 1 ∈ ∂D, v 1 points outwards the domain D and Equilibrium trajectories in the refraction billiard.Left: concatenation of non-homothetic inner and outer arc that refract one in the other; the existence of this kind of trajectory will be proved analytically, in the case of a circular domain, in Proposition 3.5.Right: examples of homothetic equilibrium trajectories.Figures taken from [28,7].final conditions (p 1 , v 1 )can be again parametrised through a pair of one dimensional quantities (ξ 1 , α 1 ), and can be used as initial conditions for a new outer arc: the above machinery can be then iterated to obtain a new concatenation.The map is called first return map, and can be used to describe the dynamics of our billiard in the phase space, parametrised by the variables (ξ, α), every time a complete concatenation of outer and inner arc is performed.At the moment, we don't make any assumption D, except for the smoothness of its boundary; on the other hand, as ordinary in billiards theory, the dynamical properties of the system (good definition of F, existence of equilibrium/periodic orbits, stability of the latters, integrable rather than chaotic behaviour) depend crucially on the geometric features of ∂D.Sections 3.2, 3.3 and 3.4 aim to describe, from different points of view, such complex interdependence between the geometry of D and the dynamics of our billiard.

3.2.
Equilibrium trajectories, stability and bifurcations.Whenever a new dynamical system is taken into account, it is quite natural to start its analysis by searching for its equilibrium trajectories, as well as investigate their stability, using the tools of nonlinear analysis.In the formalism of the first return map, equilibrium trajectories of the two dimensional system correspond to fixed points for F (see for example Figure 3.3).In the case of the refraction billiard, there is a particular class of equilibrium trajectories, called homothetic, whose existence is ensured provided very simple conditions on the boundary are verified: such trajectories result to be of paramount importance for the analysis of the model in many different circumstances.Let us suppose to have ξ ∈ I such that: (1) the position vector for the origin to γ( ξ) is orthogonal to ∂D, namely, γ( ξ)⊥ γ( ξ); (2) the segment from the origin to γ( ξ) intersects ∂D only once.
If this happens, it is easy to show that the straight half-line from the origin in the direction of γ( ξ) is invariant under both inner and outer dynamics (note that, in the case of the inner dynamics, a Levi-Civita approach to regularise the collision at the origin has been employed, cfr.[51,27]), and it is not deflected by Snell's law (see Figure 3.3, right).Along the direction defined by γ( ξ) it is then possible to construct an equilibrium trajectory, called homothetic, which corresponds to the homothetic fixed point ( ξ, 0) of the first return map F. We highlight that, although a bouncing after the collision with the central mass might seem odd from a physical point of view, the analytic continuation of inner homothetic arcs after the collision allows to study in details the local dynamics around the singularity, giving a clear portrait of orbits which are close to the collision, and then physically relevant.It is easy to observe that condition (1) in (3.3) is equivalent to require that ξ is a critical point for the function |γ(•)|, while condition (2) can be described as a star-convexity property of the domain D with respect to the direction of γ( ξ).In the following, we will refer to parameters ξ as in (3.3) as central configurations.
The first return map F is clearly infinitely-many well defined for any point ( ξ, 0), with ξ central configuration; actually, it is possible to prove that the good definition of F holds for sure locally around ( ξ, 0).
Proposition 3.1.Let us suppose that γ is at least C 1 , and let ξ be a central configuration.
Then, there exist two positive constants δ, ϵ > 0 such that the first return map F : (ξ 0 , α 0 ) → (ξ 1 , α 1 ) is well defined and differentiable in The proof of Proposition 3.1 relies essentially in showing the existence and uniqueness of the inner and outer arcs for any initial condition sufficiently close to ( ξ, 0) in the phase space and, in the case of the inner dynamics, the transversality of such arcs; the details are given in [27,Sections 1.3,1.4].We stress that, although well defined, at this stage we do not have an explicit expression for F, even in a neighborhood of the homothetic fixed point, since such expression depends on γ.Since F is locally differentiable around homothetic fixed points, it is natural to continue their analysis by investigating their linear stability, asking whether it depends on the geometrical properties of ∂D around γ( ξ) as well as on the physical parameters E, h, ω, µ.Such analysis can be carried on by considering the Jacobian matrix of F in such points, given by Although the explicit expression of F is not known, using the implicit function theorem and knowing the analytic expression of the homothetic solutions it is possible to obtain a closed formula for its Jacobian: such expression is given by where k( ξ) denotes the curvature of γ in ξ (see [31]).The analytic expression of DF( ξ, 0) is quite complicated, but it is easy to notice as it depends on both the geometric properties of γ up to second order and the physical parameters E, ω, h, µ.Let us note that, when D is a circle centered at the origin with radius R = |γ( ξ)|, its curvature is always equal to 1/R: in such case, the terms ϵ E and ϵ I disappear, and the Jacobian reduces simply to the identity matrix.This fact is not surprising, as it is consistent with the fact that the circular refraction billiard represents an integrable and highly degenerate case (see also Section 3.3): as a matter of fact, on a circular domain every radial initial condition (i.e.any point (ξ, 0) in the (ξ, α)−plane) defines a homothetic equilibrium trajectory, so that the homothetic fixed points are not isolated anymore, and form instead a straight line of fixed points in correspondence of α = 0.As for the general case, one can notice that the curvature of γ plays a role only in the ϵ E/I terms: for this reason, such terms can be considered corrections induced by the geometry of γ with respect to the circular case.
The linear stability of ( ξ, 0) as fixed point of F can be inferred by the eigenvalues of DF( ξ, 0) (see [42]), which we call λ 1 and λ 2 .Since F is area preserving, it holds that det DF( ξ, 0) = λ 1 λ 2 = 1, and, whenever λ 1 = λ −1 2 ∈ R, the fixed point is an unstable saddle; on the contrary, if the two eigenvalues are complex and conjugated, the homothetic fixed point is a stable center.To distinguish between the two cases, in the two-dimensional case a simple criterion can be adopted: denoted by ∆ the discriminant of the characteristic polynomial associated to DF( ξ, 0), one has that The case ∆ = 0, corresponding to λ 1 = λ 2 = 1, is highly degenerate: this is what happens for example in the circular case, and, in general, nothing can be said on the linear stability.Starting from Eq. (3.4), it is possible to give an explicit formula for the discriminant ∆, which   1) for the same parameters' value.Figures taken from [27]. is given by The sign of ∆ can be investigated numerically whenever one has an explicit expression for the curve γ: in the following, we propose a thorough illustration of the elliptic case, which, in the framework of the mathematical billiards, represents a case study of great importance (see for example [67,44]).
The elliptic case: analytical and numerical results.To give a practical example on how Eq.(3.5) can be used to give exact information on the stability of the homothetic equilibrium trajectories (and, in some cases, on the overall properties of the first return map), let us suppose that γ describes an ellipse with center in the origin, semimajor axis equal to 1 and eccentricity e, that is, In this case the only four homothetic trajectories are in correspondence of ξ(0) = 0, ξ(1) = π/2, ξ(2) = π and ξ(3) = 3π/2, and they are pairwise symmetric.For any of the corresponding homothetic points, it is then possible to compute DF( ξ(i) , 0), i = 0, . . ., 3, and, consequently, the discriminants ∆ (0) , . . ., ∆ (3) .The explicit expressions of these quantities, as well as rigorous asymptotic analysis, is provided in [27, Section 1.6]; here, we limit ourselves to an example, which is of particular significance to show the consistency between the analytical tools and the numerical results.Let us take for example the numerically computed values of ∆ (0) and ∆ (1)  displayed in Figure 3.4, left, where we fixed E = 2.5, ω = √ 2, µ = 2, e = 0.1, and the inner energy h varies in [0, 150].It is clear that, while the homothetic in ξ(0) = 0 is always a saddle, the stability of ξ(1) = π/2 changes when h increases: in the literature, as for example in [42], θ δ Figure 3.5.Left: example of brake two-periodic trajectory.Right: construction of the free fall map: given a direction defined by θ, it returns the angle δ between the refracted outer arc and the corresponding radial direction.Figures taken from [27].
phenomena where the dynamical properties of a map (stability of the fixed points, number of the latters, etc.) are referred to as bifurcations.We can compare Figure 3.4, left, with Figure 3.4, right, which shows the Poincaré map (that is, the representation in the (ξ, α)−plane of the iterates for F for different initial points) in a neighborhood of (π/2, 0) for different values of h close to the value of h at which ∆ (1) changes sign.One can clearly see as the fixed point, which initially is a center, changes its stability, becoming a saddle, and, for increasing values of h, leading to the formation of a new nonhomothetic, 2−periodic point.
The homothetic equilibrium trajectories analysed up to now are of great importance also for the further analysis, and in particular in Section 3.4; nevertheless, there exists another class of (two-periodic) equilibrium trajectories whose existence can be derived by purely analytic arguments.This is the case, for example, of the two periodic brake orbits composed by a pair of outer homothetic arcs connected by an inner hyperbola (see Figure 3.5, left).Such kind of trajectories can appear whenever an inner arc refracts in both sides in radial directions, and their existence can be showed analytically by means of a free fall method.The general idea under the free fall in this case is to construct a univariate function that, given the initial conditions corresponding to an homothetic outer arc, follows the generated trajectory until it exits again from D and returns the angle between the subsequent outer arc and the radial direction in the exit point (Figure 3.5, right).In this way, provided the above function (called free fall map) is well defined (and, possibly, differentiable), the search for two periodic brake orbits translates in searching for zeroes of a continuous function.The good definition of the free fall map follows from a more general geometric property of elliptic boundaries.Since this result is interesting also by itself, we write it down in the following Proposition.Proposition 3.2.[27, Proposition 6.3] Let D be an elliptic domain whose boundary is parametrised as in (3.6), with e ∈ [0, 1/ √ 2).Then, for any E, h, µ > 0, every Keplerian arc of energy E + h and mass parameter µ intersects ∂D at most in two points.
Let us remark that 1/ √ 2 ≃ 0.707: the above Proposition holds then for a wide class of ellipses, not necessarily close to a circle.Whenever e ∈ [0, 1/ √ 2) the free fall map can be proved to be well defined, the existence of brake two-periodic trajectories for suitable values of the physical parameters can be proved.Theorem 3.3.[27, Theorem 6.4] Fixed every E, ω > 0 such that ω 2 > E and any ellipse with the center at the origin, semimajor axis equal to 1 and e ∈ [0, 1/ √ 2), if µ and h are sufficiently large, then the first return map admits at least four two-periodic brake trajectories.
The results proposed until now hold in a local sense for quite generic domains and, in a more global setting, in the special case of elliptic domains.In Section 3.3 we will try to provide global results for a more general class of domains: it is the case of the close to circle ones, which will be analysed through the powerful tools brought forth by perturbation theory.

3.3.
Quasi-circular domains: a perturbative approach.It is already clear from the discussion in Section 3.2 that the shape of the domain D is of fundamental importance to infer the properties of the billiard map; in the circular case, this become evident, as the central symmetry of the domain has radical consequences on F, which will be discussed in the next paragraph.The system, in such case, results to be globally well defined, completely integrable and admits orbits of any rotation number within a certain interval.When the domain D is sufficiently close to a circle, one can ask whether some of these properties (as for example the good definition or the existence of some particular orbits) are still maintained: in this Section we aim to answer this question, taking advantage of tools coming from perturbation theory and general facts holding for area preserving maps (for a wide dissertation on the subject, see [39]).Such instruments require a definition of the analytical framework we are working in which is a bit deeper than the one described in the previous Sections, in particular involving the so-called generating function (see [39]), that now we will define in our specific case.Let us assume again that the boundary of our domain can be parametrised by a curve γ: for the sake of simplicity, we will assume that γ is 2π-periodic, and, with an abuse of notation, we still denote with γ the periodic extension of the curve, namely, γ : R /2πZ → R 2 , where R /2πZ denotes the 2π−periodic torus.Let us now take the function where L E and L I are the outer and inner Jacobi distance as defined in Eq. (3.2).By means of the implicit function theorem, the parameter ξ can be expressed as a function of ξ 0 and ξ 1 from the relation provided the following non degeneracy condition holds Recalling the variational interpretation of Snell's law, one can notice that, given two points p 0 = γ(ξ 0 ) and p 1 = γ(ξ 1 ), the generating function G returns the Jacobi length of the concatenation that connects p 0 to p 1 with an outer and inner arc, and trespasses the boundary precisely at the point p = γ( ξ) which ensures that the refraction law is satisfied by the arcs.The good definition of G, ad well as its differentiability, is not always guaranteed, and depends on D; more precisely, it is associated to the existence and uniqueness of outer and inner arcs connecting any pair of points on ∂D, and must be verified case by case.
Generating functions are commonly used when dealing with billiards (see also [66]), since the first return map, in a suitable set of canonical action-angle variables, can be implicitly expressed in terms of derivatives of G.In particular, one can define the canonical actions conjugated with the parameters ξ 0 , ξ 1 : (3.9) such quantities, which in the following will replace the angles α 0 , α 1 in the construction of a first return map, have in turn a geometrical interpretation, given by (see [28]) Eq.(3.9) translates in the fact that, whenever the initial and final point of a concatenation are known, the initial and final actions I 0 and I 1 (and, as a consequence, the angles α 0 and α 1 ) can be computed from the generating function's derivatives.Starting from this, it is possible to reconstruct the first return map in terms of the variables (ξ 0 , I 0 ) by means again of the implicit function theorem: whenever one can invert the first relation in (3.9) to obtain ξ 1 as a function of ξ 0 , and I 0 , and then define the first return map5 as The good definition of the first return map in suitable regions of R /2πZ × R is ensured whenever G E and G I are well defined and the nondegeneracy conditions (3.8) and (3.10) hold: such hypotheses will be verified case by case, possibly using different techniques.
The circular case.The final aim of the investigation presented in the current Section is to provide dynamical results holding for the billiard map induced over a domain which is quasi circular.To do it, we will adopt a perturbative point of view, taking as a unperturbed case the circle and then applying slight modifications to the boundary.To do this, a careful analysis on the map whenever the domain is circular is in order.Let us suppose that D is a circle of radius 1, whose boundary is parametrised by γ(ξ) = (cos ξ, sin ξ), ξ ∈ R /2πZ : in this case, both the potentials and the boundary are centrally symmetric, and, geometrically, this translates in the conservation of the angle α after every concatenation of outer and inner arc (see Figure 3.6, left).Taking advantage of this fact it is possible to give an explicit formulation for the first return map in this case, given by (3.11) ) (ξ 0 , I 0 ) = (ξ 1 (ξ 0 , I 0 ), I 1 (ξ 0 , I 0 )) = (ξ 0 + θ E (I 0 ) + θ I (I 0 ), I 0 ), where I c = E − ω 2 /2.Equation (3.11) has been obtained by means of analytical and geometrical reasonings, coming for classical results of Celestial Mechanics as well; the detailed computations can be found in [28].We can observe that F (c) is in the form of a shift map, that is, a map that at every iterate operates a change in the angle which depends only on I 0 by keeping constant the action.This is a direct consequence of the rotational invariance of the problem in the circular case, and shows as the action, once fixed by the initial conditions, is an integral of the motion of the system.For this reason, we can say that, in the circular case, the first return map is completely integrable, since it is a two-dimensional discrete map and has two conserved and independent quantities.The shift in the angle is a C 1 function of I 0 , and can be splitted into two terms: θ E represents th angle displacement after the outer arc, while θ I is the shift in the final angle after the Keplerian inner arc (see Figure 3.6).We can clearly say that F (c) is globally well defined in R /2πZ × (−I c , I c ); moreover, one can check that, for any I ∈ (−I c , I c ), one has that θ ′ E (I) > 0 and θ ′ I (I) < 0. The first problem we want to address in the circular case is the possible presence of periodic orbits; to do this, let us introduce the concept of rotation number (see [39]).Although the above definition holds in general, it is clear that, in the circular case, the rotation number depends only on I 0 and is simply given by the shift θ E (I 0 ) + θ I (I 0 ).Orbits with 2π−rational rotation number can be periodic, in the sense that, supposing ρ(ξ 0 , I 0 ) = 2πp/q, one has a (p, q)−periodic orbit such that (ξ q , I q ) = (ξ 0 + 2πp, I 0 ) = (ξ 0 , I 0 ); on the other hand, whenever ρ(ξ 0 , I 0 ) / ∈ Q, the corresponding orbit covers densely the invariant line {(ξ, I 0 ) | ξ ∈ R /2πZ } in the (ξ, I)−plane (see Figure 3.7).The line {(ξ, 0) | ξ ∈ R /2πZ }, although still invariant, represents an exception to the dichotomy between periodic (but moving) and dense orbits: it corresponds to all the initial conditions for the homothetic fixed points (see Section 3.2), and it is then covered in a continuum of points with rotation number equal to 0. The following Proposition resumes the results obtained in the circular case in terms of existence of (periodic and non-periodic) orbits with fixed rotation number; for the analytical expression of all the threshold values, as well as the proof, one can check [28].Proposition 3.5.[28, Propositions 4.9 & 4.10] There exists a constant C ∈ (0, π), which depends on the physical parameters E, h, µ and ω, such that for every ρ ∈ (−C, C) there exists an action value I ∈ (0, I c ) for which, for every ξ 0 ∈ R /2πZ , one has ρ(ξ 0 , ±I) = ρ.In particular, for any p, q ∈ Z such that 2πp/q ∈ (−C, C), there exist at least two (p, q)−periodic orbits starting at any ξ 0 ∈ R /2πZ .Furthermore, for an open set of values of the physical parameters (that is, an open set in the (E, ω, µ, h)−space opportunely defined), there exist at least two non-homothetic fixed points of F (c) .The orbit lie on horizontal invariant lines with I = const, and could be either periodic (the dotted lines) or cover densely the line (examples are the continuous lines on the top and bottom of the figure).A particular case is given by the invariant line I = 0, which is covered in homothetic fixed points (see Section 3.2). Figure taken from [28].
Note that the second part of Proposition 3.5 proves analytically the existence of non-homothetic fixed points in the circular case, as the one displayed in Figure 3.3, left.Quasi-circular domains.Once the analysis of the unperturbed case has been carried on, one can ask whether and under which conditions Proposition 3.5 continues to hold when we are dealing with domains which are close to be circular.The class of the quasi-circular domains we are taking into account is obtained by performing a radial deformation of the boundary depending on a generic smooth function and on a parameter ϵ, namely, where f (ξ, ϵ) is a one-valued function smooth in both variables and C ϵ > 0 is arbitrarily large.In this way, for any ϵ > −1 the parameter ξ still represents the angle between the corresponding point γ(ξ; ϵ) and the x−axis.We will refer as D ϵ to the domain whose boundary is parametrised by γ(•; ϵ).In general, when ϵ ̸ = 0 the central symmetry which characterise the circular case breaks, and we are no longer able to provide an explicit expression for the first return map on ∂D.Nonetheless, it is possible to prove an analogous of Proposition 3.5, by using a more powerful set of tools, coming from the general theory of area preserving maps.To take advantage of them, it is necessary to consider again the generating function as defined in (3.7): where again ξ can be implicitly defined in neighborhoods of points (ξ 0 , ξ 1 ) for which condition (3.8) is satisfied.Restricting the domain in the actions and assuming that the perturbing function f in (3.13) is regular enough, as well as ϵ small enough, on can prove that the first return map generated by G(•, •; ϵ) is well defined, and can be expressed as a (unknown, in principle) perturbation of (3.11).The need of restricting the action domain is related to the nondegeneracy conditions (3.8) and (3.10).In general, two-dimensional area preserving maps for which (3.14) holds are called twist maps: such property will be important for the forthcoming analysis.By continuity with respect to ϵ, whenever Proposition 3.6 holds the map F(•, •; ϵ) can be written as a perturbation of F (c) of the form where A and B are two unknown functions of class C k−2 and tend to the constant zero-function whenever ϵ → 0 in the C k−2 −norm.Now that we have defined the first return map, as well as its domain, for quasi-circular billiards, we can build up a proving algorithm, which, taking advantage of some general results for area preserving maps, will ensure the existence of orbits with fixed rotation number also in the perturbed setting.To do this, we will use two powerful results of nonlinear analysis, namely, the Poincaré-Birkhoff theorem and the Aubry-Mather one.In general (for a more rigorous explanation, see [39,28]), these theorems apply to an area-preserving twist homeomorphism on an annulus R /2πZ × [a, b], which preserves the boundaries R /2πZ × {a} and R /2πZ × {b} with rotation numbers ρ a and ρ b .As for Poincaré-Birkhoff theorem, it focuses on the periodic case, claiming that for any 2π−rational number ρ ∈ (ρ a , ρ b ) there are at least two periodic orbits whose rotation number is precisely ρ.Aubry-Mather theorem extends the result to any real number in the interval (ρ a , ρ b ), claiming the existence of at least one orbit for any of them.Before applying these results, it is necessary to construct an invariant set of F(•, •; ϵ) where the map is surely well defined, area preserving and twist and whose boundaries are invariant curves for F. To do this, we shall take advantage of the celebrated KAM theorem (see [55]), an extremely powerful result in perturbation theory that ensures the persistence of suitable invariant curves under small changes in a shift map as in (3.11).The overall proving scheme towards the final result can be resumed as follows: • we construct an invariant set K by means of KAM theorem.More precisely, we prove that for ϵ sufficiently small there exist two curves over the ξ−axis which are invariant under F(•, •; ϵ) and have irrational rotation number.The set K is precisely the region of the (ξ, I)−plane bounded by these two curves (see Figure 3.8); • with a suitable change of coordinates, we will deform K to put it in a form of an annulus R /2πZ × [a, b]: in the new coordinates, the perturbed map F satisfies the hypotheses of both Poincaré-Birkhoff and Aubry-Mather theorems; • we apply the two theorems to obtain, finally, our existence result.The most delicate part in the above scheme consists in finding a regime (namely, a region in the (ξ, I)−plane) and a sufficiently small value of ϵ for KAM theorem to be applied.As a first point, we shall choose carefully the curves which we want to preserve under perturbation: they must have a Diophantine rotation number (see [55] and [28]).We also stress that, in view of Proposition 3.6, the above procedure can be applied to any connected component of (−I c , I c ) \ J , so the multiplicity of the found orbits can change accordingly.We present now the final result of this Section without any further detail on the proof: nonetheless, we invite the interested reader to go through the application of all the proposed techniques in [28].− and θ i + as the rotation numbers for the unperturbed dynamics for I 0 = a i and I 0 = b i (for simplicity, let us assume that θ i − < θ i + ).For any i = 1, . . ., N , let us fix ρ i ± two Diophantine numbers such that θ i − < ρ i − < ρ i + < θ i + .Then there exists ε > 0 such that, for every ϵ ∈ R, |ϵ| < ε, and every ρ ∈ (ρ i − , ρ i + ) there is at least k orbits for the perturbed map with rotation number ρ, where k is the number of intervals (ρ i − , ρ i + ), for different i, in which ρ is contained.Moreover, if ρ is 2π−rational, then the orbits are 2k and they are periodic.
Although a little bit technical in its set up, the core of the theorem is that, restricting suitably the set of rotation numbers and taking domains which are sufficiently regular and close to a circle, the existence of orbits with fixed ρ, including periodic ones, is guaranteed.The existence of a wide variety of periodic orbits for a non-circular case is a highly nontrivial results, and can be interpreted as a strong hint of the presence of a complex dynamics, which might be chaotic: this is precisely the topic of Section 3.4, where we will search for conditions on the domain's shape under which it is possible to prove analytically the chaoticity of the model.

3.4.
The onset of chaos in galactic refraction billiards.The study of the dynamics of galactic refraction billiards, carried on starting from the existence and stability of equilibrium trajectories in Section 3.2 and continued with the analysis of the quasi-integrable setting induced on the model by the choice of a quasi-circular boundary, finds in this last Section its conclusion.Here, the problem addressed is the possible chaoticity of our model, and, in particular, the existence of simple, geometrical conditions on D that ensure that the system satisfies a mathematically rigorous definition of chaos.Such question, which in some sense was the first reason to motivate us in investigating the galactic refraction model, comes quite natural while observing some of the simulations provided in [27] (see also Figure 3.9).At least in the particular case of an ellipse having its center at the origin, it is evident how increasing values of the inner energy h lead to the appearance of diffusive orbits around the homothetic points on the ξ−axis, providing a clear evidence of chaotic behaviour, in a subset of the phase space.
Of course, numerical simulations as the ones presented in [27] are not enough to prove the actual chaoticity of the model, as numerical instabilities can get in the way and lead to possibly non-accurate deductions.It is then important to go further with a rigorous proof, to show analytically how chaotic phenomena can occur under precise hypotheses.A first problem one has to deal with is which definition of chaos it is convenient to take into In the first case, a circle whose center is at the origin admits infinitely-many central configurations, but none of them are nondegenerate.On the other hand, an ellipse with one of the foci at the origin admits two nondegenerate central configurations, but they are antipodal.
account: roughly speaking, a dynamical system is considered chaotic whenever it is sensitive to changes in the initial conditions.More precisely, when, moving from a given point, the trajectories' behaviour becomes unpredictable, potentially covering the whole phase space.From a mathematical point of view, there are many different ways to explicitly describe such concept (see for example [41]): in the current paper, we will use Devaney's definition of topological chaos, presented in full details in [30].
The main result regarding chaotic behaviour in a galactic refraction billiard, presented in [5], resumes in finding a geometrical condition on the boundary, called admissibility, which ensures the existence of a topologically chaotic subsystem for high enough inner energies.
Definition 3.8.Let us take a domain D containing the origin, and, taking γ : I → R 2 as a parametrisation of ∂D, suppose that it is at least C 2 .The domain D is termed admissible if there are at least two nondegenerate and non-antipodal central configurations, namely, if there exist ξ 1 , ξ 2 ∈ I as in Eq. (3.3) which are strict maxima or minima for the function ∥γ(•)∥ and such that γ(ξ 1 ) and γ(ξ 2 ) are not collinear with the origin (see Figure 3.10).
Theorem 3.9.[5, Theorem 4.7] Let D be an admissible domain as in Definition 3.8.Then, if the inner energy's parameter h is high enough, the galactic refraction billiard is chaotic, in the sense that it admits a topologically chaotic subsystem.
While in [7] the proof of Theorem 3.9 is explained in details, here we limit ourselves in pointing out that such chaotic subsystem is obtained by conjugation with the Bernoulli shift, achieved by constructing a suitable symbolic dynamics (see [30]), which is well defined in a subset of the initial conditions and whenever the inner energy is high enough with respect to the outer one.For the sake of completeness, let us briefly resume the main ideas behind such construction.In general, we say that a dynamical system F : A → A, A being the set of initial conditions, admits a symbolic dynamics whenever there exists a surjective and continuous projection map Π : A → {1, . . ., n} Z , n ≥ 1, such that the diagram A A {1, . . ., n} Z {1, . . ., n} Z F Π Π σ commutes, where the function σ is the Bernoulli shift: given a sequence in {1, . . ., n} Z , σ moves any of its elements on the right 6 .In practice, constructing a symbolic dynamics corresponds to encoding the forward and backward orbits F into sequences of symbols, whose geometrical and physical meaning depends by the model itself.This fact is particularly important when Π is bijective: in this case, we say that our system is topologically conjugated with σ, and the one-to-one correspondence between F -orbits and bi-infinite sequences of symbols allows to conclude that F is chaotic.
In our case, our aim is to construct a topological conjugacy between the first return map F : (ξ 0 , α 0 ) → (ξ 1 , α 1 ), possibly restricted to a subsystem, and the Bernoulli shift on a suitable set of sequences.The map Π keeps tracks of the points where subsequent concatenations of outer and inner arcs, and in particular the ones that remain close to the homothetics, intersect the boundary.With reference to Figure 3.11, we proceed as follows: given an admissible domain, suppose that it admits m central configurations ξ 1 , . . ., ξ m , where m ≥ 2 by definition.
We then construct U 1 , . . ., U m ⊂ ∂D suitable neighborhoods of γ(ξ 1 ), . . ., γ(ξ m ), and say that a trajectory starting with initial conditions (ξ 0 , α 0 ) realizes the word s ∈ {1, . . ., m} Z if and only if it crosses ∂D only in ∪ m i=1 U i and in the order prescribed by s itself: in such case, we define the projection Π(ξ 0 , α 0 ) = s.Putting some restrictions on the neighbourhoods U i as well as on the bi-infinite sequences in {1, . . ., m} Z , one can prove that: • there exists a set X of initial conditions (ξ 0 , α 0 ) for which Π is well defined, that is, the backwards and forwards trajectories with initial conditions (ξ 0 , α 0 ) ∈ X cross ∂D close to a central configuration (it is said that such trajectories shadow the homothetic ones); • under the hypotheses of nondegeneracy and non-antipodality of central configurations and for h sufficiently large, Π is bijective and continuous on X, and hence F | X is topologically conjugated to Bernoulli shift, thus chaotic.To prove the surjectivity of the map, we shall use of Poincaré-Miranda theorem (see [53]), a topological fixed point result, and of the variational characterisation of refracted arcs.In [7], a thorough analysis on the admissibility condition, including results that still hold when some of the requirements are weakened, is carried on; furthermore, conditions on the bi-infinite sequences which lead to a non-collisional dynamics (namely, whose trajectories do not collide with the central mass) are detected.
Remark 3.10.Although [7] is focused on the case of a refractive galactic billiard, it can be proved that a result analogous to Theorem 3.9, with the same admissibility conditions, holds in the class of Kepler reflective billiards .In this case, the chaoticity result represent an interesting complement of the result obtained by Takeuchi and Zhao in [67].Here, as a consequence of a more general result, they prove that a Kepler reflective billiard whose boundary is a focused ellipse (i.e. an ellipse with one of the foci put at the origin, see Figure 3.10, Left) induces an integrable dynamics.
On one hand, Takeuchi and Zhao's result is completely coherent with Theorem 3.9, since a focused ellipse is not an admissible domain according to Definition 3.8 (it indeed has two nondegenerate central configurations, parallel to the x−axis, but they are antipodal); on the other hand, the presence of a transition between integrable and chaotic systems by simply translating the boundary is a nontrivial fact that is worthy of a further investigation.In [6], we start by proving, analytically, that elliptic (galactic refraction and Kepler reflection) billiards are almost always chaotic, provided the inner energy is sufficiently large.More precisely, fixed a non-circular elliptic shape with any eccentricities, it can be put almost everywhere with respect to the origin (keeping 0 in its interior) to have an admissible domain.
We conclude this Section by highlighting how admissibility is only a sufficient condition ensuring chaos for large inner energies: it is then natural to ask what happens when such condition is violated but there are no evidences of integrability.In such case, one may ask themselves whether the limit of our proof is only constructive, in the sense that, taking a different projection map, some of the admissibility conditions can be weakened while still obtaining chaos.This will be the subject of a forthcoming paper.

Conclusions and further perspectives
This paper aimed to show the potential of analytical and semi-analytical techniques when applied to the investigation of dynamical systems coming from Celestial Mechanics, by presenting two examples, the first one of which was on a planetary scale and the second one on a galactic one.
The first problem addressed was the search of stability estimates for the motion of a small body moving around the Earth in the so-called geolunisolar model, in the sense that stability times up to which the variation of the orbital elements of a given geocentric orbit have been computed.The results have been obtained by means of two different techniques: the first, based on normal form theory, provided stability times of the order of 10 4 years for the quantity 1 − e 2 (1 − cos i), e and i being respectively the eccentricity and the inclination of the orbiting body, which hold for quasi-circular and quasi-equatorial initial trajectories.The second method, based on the celebrated Nekhoroshev theorem, finds its application on a wider set of initial conditions in terms of eccentricities and inclinations; as for the semimajor axis, the results in this case hold for a strip in MEO from 11000 to 19000 km, with stability times that, starting from being very long (10 5 years) for low distances, tend to decrease with the altitude.The reason for such worsening in the estimates can be found in the convergence of the constructive algorithm used, and could be potentially improved.On the other hand, another fact which worsen our estimates are given by the presence of lunisolar resonances: in such case, a certain improvement could be achieved by considering Nekhoroshev theorem in its complete form instead of the nonresonant one, and proceeding with a careful analysis of our secular geolunisolar Hamiltonian whenever a resonance occurs.We stress that the approach presented in Section 2 could be potentially used to provide stability estimates in any gravitational system where small masses orbit around a central extensive body, of which the approximated shape (and, as a consequence, the corresponding gravitational potential) is known, even in presence of third bodies which act as a perturbation.
As for the second model considered, it has been studied with analytical techniques coming from billiards' and, more generally, classical dynamical systems' theory.In this case, many results regarding the presence of equilibrium, periodic and quasi-periodic orbits are provided, along with evidences that, under suitable conditions, the central mass acts as a scatterer, deflecting our particle in a chaotic, an somehow unpredictable, way.In general, refraction billiards where two differential potentials are coupled through a refraction interface can be used as a simplified model to study any complex dynamical system whose behaviour presents different regimes, for example depending on the particle's position; in such general case, the same formalism, as well as the same techniques, presented in Section 3 can be used.

Figure 1 . 1 .
Figure 1.1.Examples of orbits of refraction galactic billiards.The orbit goes inside and outside the domain, being deflected at every passage through the interface.Left: three-periodic trajectory.Right: quasi-periodic trajectory (figure taken from[28]).

Figure 2 . 1 .
Figure 2.1.Stability times computed for different values of semimajor axis, eccentricity and inclination using the nonresonant version of Nekhoroshev theorem.The color scale refers to the computed stability times (in years), while the white region correspond to the values of (e, i) where Theorem 2.1 can not be applied with the present algorithm.The red lines are in correspondence of the inclinations of the secular geolunisolar resonances.Data taken from [15].

3. 1 .
The model: analytical set-up and motivations.Let us start the analytical description of our galactic refraction billiard by taking a smooth open domain D ⊂ R 2 , containing the origin, and considering the potential

1 pFigure 3 . 1 .
Figure 3.1.Left: Snell's refraction law.The angles α E and α I are the angles respectively of the outer and the inner arc with respect to the normal direction to ∂D in z.The two angles are connected by the relation (3.1).Right: concatenations from p 0 to p 1 with an outer and inner arc, for different positions of the transition point p.The left figure is taken from [7].

Figure 3 . 6 .
Figure 3.6.Left: conservation of the angle α in the circular case.Right: Outer and inner shift in the circular case.By virtue of the central symmetry of the overall system, they depend only on I 0 and not on ξ 0 .

Figure 3 . 7 .
Figure 3.7.Orbits of the refraction galactic billiard in the circular case, taking into account the physical parameters' values E = 7, ω 2 = 3, h = 2 and µ = 15.The orbit lie on horizontal invariant lines with I = const, and could be either periodic (the dotted lines) or cover densely the line (examples are the continuous lines on the top and bottom of the figure).A particular case is given by the invariant line I = 0, which is covered in homothetic fixed points (see Section 3.2). Figure taken from[28].

Figure 3 . 8 .
Figure 3.8.Construction of the invariant set K for the perturbed first return map F(•, •; ϵ) in the (ξ, I)-plane.The curved solid lines are the invariant curves for the perturbed dynamics, whose existence is ensured by KAM theorem and with rotation number ρ a and ρ b .The dashed straight lines are the corresponding orbits in the unperturbed dynamics.

Table 1 .
Stability time (in years) for the quantity I