Camille Saint-Martin

Adrien Morel

Ludovic Charleux

Emile Roux

Aya Benhemou

Quentin Demouron

Adrien Badel

Escape from the potential well of bistable vibration energy harvesters using buckling level modifications

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The availability of vibrations in our environment and the development of wireless sensor networks motivates researchers to develop energy harvesting systems to supply low-power electronic systems [START_REF] Prajwal | A review on vibration energy harvesting technologies: analysis and technologies[END_REF]. NVEH exhibit broader bandwidths compared to their linear counterpart, making them particularly relevant for harvesting energy from broad and time-varying vibration spectrum [START_REF] Sebald | Simulation of a Duffing oscillator for broadband piezoelectric energy harvesting[END_REF]. However, NVEH exhibit complex dynamics i.e., there exists multiple orbits of different power for a given vibration frequency. Some orbits are of low power, which may lead to poor harvesting performances. Therefore, to improve the performance of NVEH, researchers worked on the development of orbit jump strategies to shift from low-power orbits to higher power orbits [START_REF] Huguet | Orbit jump in bistable energy harvesters through buckling level modification[END_REF]. In this study we use a bistable NVEH based on a buckled beam to which we apply an orbit jump strategy based on the modification of the buckling level.

Results and discussion

Figure 1(a) shows a schematic representation of the Duffing-type piezoelectric vibration energy harvester with two stable equilibrium positions (EP) for x = ±x w considered from [START_REF] Huguet | Orbit jump in bistable energy harvesters through buckling level modification[END_REF]. This NVEH is composed of a mechanical resonator and an electrical extraction circuit. The NVEH is excited with a sinusoidal acceleration of amplitude A and frequency f d . The harvested power is the one dissipated in the extraction circuit, which is simply a resistor R. This orbit jump consists in dynamically modifying the buckling level by adjusting the EP x w of the NVEH. We supposed that the system starts on a low-power orbit. The initial value of the EP x w (0) is x w , thereafter we modify the EP with a factor k 0 w at t 0 : x w (t 0 ) = k 0 w x w (0). The new value of the EP is maintained for a certain duration ∆t. Then, from the instant t 0 + ∆t the EP is further modified with a factor k 1 w such as: x w (t 0 + ∆t) = k 1 w x w (0). Figure 1(a) shows the evolution of the EP. Therefore, the investigated orbit jump strategy depends on four influence parameters (t 0 , ∆t, k 0 w , k 1 w ). We used an evolutionary algorithm in order to find optimum parameters combination. The criterion to maximize is the total harvested energy on a duration of 100 T d that takes into account the invested energy to modify the buckling level. Figure 1(b) shows an example of a successful application of such an orbit jump strategy for f d = 50 Hz. w , k 1 w ) = (2 × 10 -3 , 7 × 10 -3 , 5.30, 3.50). The system starts on a low power orbit (blue) and converges to a high power orbit (red).

Figure 1 :

 1 Figure 1: (a) NVEH scheme. (b) Trajectory in the dimensionless phase plane for f d = 50 Hz, A = 2.5 m/s 2 with (t0, ∆t, k 0w , k 1 w ) = (2 × 10 -3 , 7 × 10 -3 , 5.30, 3.50). The system starts on a low power orbit (blue) and converges to a high power orbit (red).


Frédérique Charles
email: frederique.charles@sorbonne-universite.fr

Annamaria Massimini
email: annamaria.massimini@tuwien.ac.at

Francesco Salvarani
email: francesco.salvarani@unipv.it

Mathematical and numerical study of a kinetic model describing the evolution of planetary rings

In this paper, we study a kinetic model describing the evolution of planetary dust under the action of a planet and its satellites. In particular, we focus our attention on the formation of planetary rings and the role of shepherd moons. The equation describing the considered physical phenomenon is of Vlasov type, posed in an evolutionary domain in time. We first study some theoretical properties of the model. Then, we describe a numerical method, suitable for the study of kinetic equations in evolutionary domains with possibly complicated geometries. Finally, we show and comment some simulations. In particular, our numerical simulations show that shepherd moons play a key role in the formation and maintenance of divisions between rings.

Introduction

A planetary ring is a complex system composed of dust and other small particles that orbit a planet forming a flat disk. The first observed planetary rings were those of Saturn and were recognized as rings by Christiaan Huygens in 1655 [START_REF] Huygens | Systema Saturnium, sive de causis mirandorum Saturni phaenomenon, et comite ejus planeta novo[END_REF]. However, ring systems are not a feature unique to Saturn. In fact, the other three giant planets of the Solar System are also surrounded by a system of rings, and astrophysicists infer that many exoplanets may also have ring systems [START_REF] Kenworthy | Mass and period limits on the ringed companion transiting the young star J1407[END_REF].

In 1676, Giovanni Cassini discovered a gap between the rings of Saturn, now called the Cassini Division [START_REF] Domenico | Observations nouvelles touchant le globe & l'anneau de Saturne[END_REF]. Thanks to the Voyager probes, it was discovered that the structure of the rings is very complex. Sometimes, within gaps in the rings are moons, called shepherd moons. The gravity of the shepherd moons serves to maintain a well-defined edge of the ring. Material approaching the orbit of the shepherd moon can be deflected back into the ring body, ejected from the system, or fall onto the moon itself.

The formation of rings and gaps between rings has logically attracted the attention of astronomers and astrophysicists (see, for example, [START_REF] Esposito | Planetary rings: a post-equinox view[END_REF] and the references therein).

As a natural consequence, this problem has also been studied by applied mathematicians. In this context, multiple approaches have been proposed. The first one consisted in 1 considering the problem from the microscopic point of view, by using the classical methods of celestial mechanics (see, for example, [START_REF] Felino | On periodic perturbations of uniform motion of Maxwell's planetary ring[END_REF][START_REF] Benet | Phase-space structure for narrow planetary rings[END_REF]). Always at microscopic level, the literature has also been interested in the study of models in which the particles constituting the ring suffered not only the effects of the gravitational field, but could also collide inelastically [START_REF] Kawai | An inelastic collision model for the evolution of "planetary rings[END_REF].

The second viewpoint is the mesoscopic (or kinetic) approach. It consists in describing the system through particle densities in the phase-space of the system. This approach goes back to James Clerk Maxwell [START_REF] Garber | Subjects great and small: Maxwell on Saturn's rings and kinetic theory[END_REF][START_REF] Felino | On periodic perturbations of uniform motion of Maxwell's planetary ring[END_REF], one of the founders of kinetic theory. A recent application of this description to planetary rings can be found in [START_REF] Majorana | Approximate explicit stationary solutions to a Vlasov equation for planetary rings[END_REF], which studies approximate explicit stationary solutions in the context of a non-collisional model of Vlasov type.

The third strategy consists in treating planetary rings by using an hydrodynamic description (see, for example, [START_REF] Spahn | Hydrodynamic description of planetary rings[END_REF]).

We moreover mention the possibility of using stochastic tools, such as in the case of the Burgers-Zeldovich model [START_REF] Neate | A stochastic Burgers-Zeldovich model for the formation of planetary ring systems and the satellites of Jupiter and Saturn[END_REF].

Independently on the scale used for describing the problem, many authors have focused themselves on a particular phenomenon occurring inside planetary rings. For instance, the literature reports studies on fragmentation and coagulation [START_REF] Esposito | A predator-prey model for moon-triggered clumping in Saturn's rings[END_REF] or on the fractal structure of rings [START_REF] Malyarenko | Fractal planetary rings: energy inequalities and random field model[END_REF].

In our paper, we focus on a particular phenomenon, namely the formation of a gap between rings caused by the effect of a shepherd moon. We use the kinetic approach and, working on a short time scale, describe the phenomenon with a non-collisional Vlasov equation in a time-dependent domain. The precise assumptions that justify the mathematical structure of the model are described in detail in the next section. At this point, we recall that the mathematical study of kinetic equations in evolutionary domains is still in its early stages (see, for example, the pioneering paper [START_REF] Bardos | Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d'approximation; application à l'équation de transport[END_REF] and [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture[END_REF][START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF][START_REF] De | GPU-accelerated numerical simulations of the Knudsen gas on time-dependent domains[END_REF][START_REF] Salvarani | On the linear Boltzmann equation in evolutionary domains with an absorbing boundary[END_REF]).

Planetary rings are essentially two-dimensional objects, because of the properties of symmetry of the problem. In our article, we hence work in the phase-space R 2 × R 2 and suppose that the spatial domain of the solution of the Vlasov equation is given by the complement of the domain occupied by a planet and by its moons.

We first study and provide the theoretical framework of the problem. We then address the problem of its numerical simulation. Although the simplifications with respect to the full dynamics make the problem linear, the existence of a moving domain is a source of difficulties that imposes a nontrivial study of the problem.

We have chosen to use a particle method for its numerical implementation. It is clear that the study of kinetic equations in an evolutionary domain is not an easy task because of the high dimensionality of the problem, whose effects need to be taken into account, and to the difficult treatment of possibly complicated geometries. In our approach, we discretize the unknowns by mean of a collection of weighted smooth shape functions, which evolve in time by following the dynamics of the problem, and then we handle the possible overlapping between ring particles and the moons (or the planet). This has as a consequence the elimination of the ring particle from the domain. The numerical results show the versatility of the numerical method. In particular, the effect of a shepherd moon on ring gap formation is clearly identified.

The structure of the article is the following. We start by discussing in Section 2 the different physical phenomena taken into account or neglected, before introducing the model. In Section 3, we prove an existence and uniqueness result for this model, by applying the method of characteristics in a moving domain. Section 4 describes the numerical strategy and Section 5 presents the numerical results for two different meaningful scenarios.

The mathematical model

We consider a planet surrounded by a fixed number of moons and a dusty cloud. The goal of our model is to describe the dynamics of ring formation. The physical processes involved in the evolution of planetary rings are multiple, and taking into account all of them would lead to considerable complications in our study. However, the effects of various physical phenomena may be more or less important with respect to the considered spatial and temporal scales. For our purposes, we hence study a simplified version of the problem, which retains, however, the main physical effects governing the creation of ringlets (in particular, the interaction between particles and shepherd satellites).

The main gravitational effects

The first step is the analysis of the orders of magnitude of the various gravitational forces acting on the system. In order to be consistent with a realistic situation, we consider, for all the four giant planets of the Solar system, the Sun-planet-moons-rings subsystem and analyze the problem in detail at the quantitative level.

The negligible role of the Sun

We first observe that the main effect on the particles is the gravitational attraction of the hosting planet and that the gravitational attraction of the Sun can be neglected in first approximation. This emerges from comparing the gravitational forces exerted on the annular particles by their host planet and the Sun. They can be directly computed from the data available for the Solar system. See the details below.

Let us consider a giant planet having mass m P and a dust particle in its ring system, with mass m. Let M be the mass of the Sun. We denote with F S and F P the gravitational forces exerted respectively by the Sun and the planet on the ring particle. The general formula for computing F P and F S on a ring particle with mass m are:

F S = GmM d 2 S , F P = Gmm P d 2 P ,
where G is the gravitational constant, d S is the average distance between the mass centre of the Sun and the particle and d P is the average distance between the mass centre of the planet and the particle. Table 1 shows some data for the planetary rings of the Solar system. Because of the complex structure of the considered ring systems, we chose for d P either the radius of one of the main rings, or an average radius.

The value of the Sun's mass is M = 2.0 × 10 30 kg. The mass m of a dust particle usually ranges from 10 -3 kg to 10 6 kg [START_REF] Brilliantov | Size distribution of particles in Saturn's rings from aggregation and fragmentation[END_REF]. However, this value has no influence in the comparison of the orders of magnitude developed in this Section.

We can deduce from these data that the gravitational force exerted by the Sun on a ring particle is negligible when compared to the gravitational force exerted by the planet on the same ring particle. Indeed, we can compute and compare the attractive forces F P (in N) between a particle of mass m and a given planet:

1. for a particle of mass m that orbits around Saturn

F P = G m (5.7 × 10 26 kg) (1.2 × 10 8 m) 2 ≈ G m × (4.0 × 10 10 kg • m -2 ),
2. for a particle of mass m that orbits around Uranus

F P = G m (8.7 × 10 25 kg) (3.8 × 10 7 m) 2 ≈ G m × (6.0 × 10 10 kg • m -2 ),
3. for a particle of mass m that orbits around Neptune

F P = G m (1.0 × 10 26 kg) (6.7 × 10 7 m) 2 ≈ G m × (2.2 × 10 10 kg • m -2 ),
4. for a particle of mass m that orbits around Jupiter

F P = G m (1.9 × 10 27 kg) (1.2 × 10 8 m) 2 ≈ G m × (1.3 × 10 11 kg • m -2 ).
On the other hand, the gravitational force F S exerted by the Sun on a particle orbiting Jupiter (which is the closest planet to the Sun having a ring system) is

F S = G m (2.0 × 10 30 kg) (7.8 × 10 11 m) 2 ≈ G m × (3.2 × 10 6 kg • m -2 ) .
In the case of the other planets considered, the gravitational force between a ring particle and the Sun is even smaller. Therefore, ring particles are mainly affected by gravitational attraction of the planet and it is reasonable to neglect the effect of the Sun.

About the gravitational and contact interaction between particles

Secondly, we neglect the mutual attraction between ring particles and assume that collisions between them are unfrequent: we hence consider only the interactions between ring particles and planet and the interactions between ring particles and moons. This assumption is more delicate. Indeed, it is clear that aggregation and fragmentation phenomena play an important role in the ring dynamics [START_REF] Esposito | Planetary rings: a post-equinox view[END_REF]. However, the importance of this phenomenon has to be appreciated with respect to the time scale of our analysis and the ring density. We consider here the situation when the collision frequency is smaller than the orbital frequency. The particle collision frequency ω c in a given ring can be computed as the product of the optical depth τ , related to the density of the ring, and the angular speed Ω of the latter, i.e. ω c = τ Ω (as explained in [START_REF] Goldreich | The dynamics of planetary rings[END_REF] and [START_REF] Frank | Structure, stability and evolution of saturn's rings[END_REF]). The first consequence that can clearly be deduced from this formula is that therefore not all rings behave in the same way in terms of collisions. For optical depth around 10 -7 , as in the most tenuous Saturn's rings, ω c ≈ 10 -11 s -1 , which means, almost 1 collision between dust particles per century. Whereas for τ = 1, as in the dense A Ring, particles collide with a frequency of 10 per day. Since the orbital period of the major Saturn's moons is of the order of magnitude of days, such as Mimas (23 hours), Thetys (45.6 hours), and the time horizon of our model is of few satellite's orbits, the system can be considered as sufficiently rarefied to make such collisions unlikely [START_REF] Griv | On the stability of Saturn's rings: a quasi-linear kinetic theory[END_REF].

About the mutual gravitational attraction of N moons

In our problem, we consider N moons and study their role during the formation of planetary rings, in the short period of time following the breakup of a satellite. While, in general, the problem is difficult to be solved (it is, in fact, the celebrated N -body problem of classical mechanics), when one of the celestial body has a mass much greater than the mass of the other bodies, it is clear that the reciprocal influence of the small-mass bodies between themselves is important only in large time. When dealing with reasonably small time intervals, the gravitational attraction between the moons can be neglected and we can therefore assume that the lunar orbits are always distinct and never overlap. These assumptions are justifiable because we work on small time scales, in which the main effect is the gravitational attraction between the planet and its moons.

The idea is to compare the gravitational force exerted by the host planet on one of its moons with that which all the remaining moons exert on the fixed moon. For simplicity, we have chosen to show the calculations for Saturn, as its moons are the most massive and numerous compared to the other ringed planets in the Solar System. The moon we consider is one of the most massive, namely Thetys. First, we calculate the Saturn-Thetys gravitational force F ST and then F DT , the one between Thethys and Dione, another massive moon of Saturn. Obviously, to be precise, we should compute the gravitational force between Thetys and all the other moons of Saturn, but it is easy to see from the data that the order of magnitude of the latter is the same as F DT . This is why we omit the calculation. The reciprocal gravitational attraction force F ST between Saturn and Tethys is

F ST = G (5.7 × 10 26 kg) × (6.2 × 10 20 kg) (2.9 × 10 8 m) 2 = 2.8 × 10 20 N,
whereas the reciprocal gravitational attraction F DT , at Dione-Thetys' minimum distance, is

F DT = G (1.1 × 10 21 kg) × (6.2 × 10 20 kg) (8.2 × 10 7 m) 2 = 6.7 × 10 15 N.
One can see that F ST is 5 orders of magnitude greater than F DT (and therefore greater than the force exerted by all the more massive moons on Thetys). Therefore, the mutual gravitational attraction between Thetys and the other moons can be neglected, working over short (astronomical) time periods. In particular, this procedure is generalisable by taking any other moon in place of Thetys. And it is generalisable by choosing any other ring system.

Interactions between particles and the major bodies

The sizes of satellites and planets being much larger than the sizes of particles, the contact interaction between the former and the dust is the predominant phenomenon in our problem. In particular, in our model we suppose that the high mass of satellites and planet causes the total absorption of particles interacting with them.

Most of the intricate ring structures owe their existence to the gravitational effect of moons, without which the rings would be flat and featureless; without moons there would probably be no rings at all because thin disks of small particles gradually would spread and disperse [START_REF] Esposito | Planetary rings: a post-equinox view[END_REF]. Resonance is a key effect in maintaining ring gaps [START_REF] Esposito | Planetary rings: a post-equinox view[END_REF]. However, during the evolution of a cloud of debris for a short period of time, we only consider the collision between satellites and particles, neglecting resonance.

The gravitational force exerted by the planet on its annular system erodes and sculpts the rings, whose particles continually rain into the planet's atmosphere. This flux of annular grains disintegrate the grains themselves, reducing the lifetime of planetary rings [START_REF] Charnoz | Origin and Evolution of Saturn's Ring System[END_REF].

We hence assume in our model that the total mass of the dust cloud is much smaller than the mass of the planet and of the mass of the satellites, and that the mass of the satellites is much smaller than that of the planet. This is the case for Saturn's ring: the Cassini probe estimated that the total mass of Saturn's rings (contained mainly in rings A, B and C) is of 1.54 ± 0.49 × 10 19 kg, which represents a fraction of the mass of Mimas [START_REF] Iess | Measurement and implications of saturn's gravity field and ring mass[END_REF]. For the other annular systems of the Solar System, whose rings are more rarefied, the argument is similar.

We will therefore suppose that the dust cloud has, in the model, a negligible gravitational effect on the system and that it is possible to neglect the mutual gravitational interactions between the satellites. Thus all the resulting complexity of the N -body problem will not be addressed in our study. Consequently, the gravitational interactions considered in this article are the action of the planet on the dust cloud and on the satellites, as well as the gravitational attraction of the satellites on the dust cloud.

Equations modelling the system planet-moons-rings

In order to write our model, we first introduce the distribution function

f : R * + × R + × Ω t × R d → R +
which describes the mass density of the dust, where m ∈ R * + is the mass variable, t ∈ R + is the time variable, x ∈ Ω t ⊂ R d is the spatial position and v ∈ R d the velocity, where d ∈ N * . In practice, because of the symmetries of the problem, the most relevant case is

d = 2.
We suppose that the origin of the reference frame is the center of mass of the planet. Let r P > 0 and r i > 0 the radii of the planet and of the i-th shepherd moon respectively (i = 1, . . . , N ), where N ⩾ 1 is the total number of shepherd moons. We introduce the sets

S t i = {x ∈ R d : |x -ξ i (t)| ⩽ r i } and P t = {x ∈ R d : |x| ⩽ r P }, where ξ i (t) ∈ R d is
the position of the center of mass of the i-th shepherd moon with respect to the origin of the reference system, i.e. the center of mass of the planet. The orbits of the satellites are assumed to be known, so ξ i (t) is a datum of the problem. Owing to the fact that we work in a short-time horizon, we moreover suppose that S t i ∩ S t j = ∅ for all i, j = 1, . . . , N and S t i ∩ P t = ∅ for all i = 1, . . . , N .

The spatial domain of definition of the problem is hence given by the following open region of R d :

Ω t = (P t ) c \ N i=1 S t i .
The boundary of Ω t is hence

Γ t := ∂Ω t = ∂P t ∪   N i=1 ∂S t i   . If M ⊂ R * + , X ⊂ Ω t and V ⊂ R d , the integral I M , X, V (t) := M × X× V f (m, t, x, v) dmdxdv
represents the number of dust particles with mass m ∈ M , position x ∈ X and velocity v ∈ V at time t.

The evolution of the dust cloud is described by the gravitational Vlasov equation

∂f ∂t + v • ∇ x f - ∇ x Φ m • ∇ v f = 0, (m, t, x, v) ∈ R * + × R + × Ω t × R d (1)
where

Φ(t, x) = Gm N i=1 m i |x -ξ i (t)| + m P |x| (2)
is the gravitational potential on a particle of mass m. It is due to the planet, which has mass m P , and to the moons, with masses m i and position of their centers of mass ξ i (t) for all i = 1, . . . , N . The planet and the shepherd moons influence the motion of the particles through the gravitational forces exerted on the dust cloud. In particular, when a dust particle collides with the planet or a satellite, it is absorbed. This effect is mathematically described by supposing that the planet and the moons are absorbing moving barriers.

Let n x be the outward normal originated in x ∈ Γ t . Then, the boundary conditions on f are the following:

   f (m, t, x, v)| x∈∂S t i , (v-v i (x))•nx<0 = 0, i = 1, . . . , N , f (m, t, x, v)| x∈∂P t , (v-v P (x))•nx<0 = 0, (3)
where v i (x) is the local velocity of the point x located at the surface of the i-th moon and v P (x) is the local velocity of the point x located at the surface of the planet. For simplicity, let us define the ingoing boundary at time t as the subset of

Γ t × R d such that Σ t -:= N i=1 (x, v) ∈ ∂S t i × R d : (v -v i (x)) • n x < 0 ∪ (x, v) ∈ ∂P t × R d : (v -v P (x)) • n x < 0 . (4)
Thanks to this definition, the boundary conditions (3) can be rewritten as

     ∂f ∂t + v • ∇ x f - ∇ x Φ m • ∇ v f = 0, (m, t, x, v) ∈ R * + × R + × Ω t × R d , f (m, t, x, v)| (x,v) ∈ Σ t -= 0. ( 5
)
The problem is supplemented with a suitable initial condition. We assume that its support is contained in Ω 0 :

f (m, 0, x, v) = f in (m, x, v) if (m, x, v) ∈ R * + × Ω 0 × R d 0 otherwise, (6)
where

f in ∈ L 1 (R * + × Ω 0 × R d ) ∩ L ∞ (R * + × Ω 0 × R d ).
3 Some mathematical properties of the model

In this section, we prove an existence and uniqueness result for the gravitational Vlasov equation for planetary rings ( 1)-( 2) with boundary and initial conditions ( 5)-( 6), exploiting the method of characteristics. To recall, the method of characteristics consists of a theoretical technique to shift the focus from the analysis of a PDE to the resolution of a system of ODEs. The first step of the method is to choose certain curves (characteristic curves) along which the starting PDE becomes a system of ODEs, via an appropriate change of variable. Secondly, the classical existence and uniqueness theorems for ODEs are invoked, with the aim of reconstructing the solution of the system along these curves.

One concludes by transforming it into the solution of the PDE, taking advantage of the change of variable previously implemented (see [START_REF] Evans | Partial differential equations[END_REF] for more explanations). In our analysis, we need to pay attention to the evolution of the spatial domain, whose changes over time are due to the motion of the planet and the shepherd satellites.

The method of characteristics for the Vlasov equation

The gravitational Vlasov equation is a linear scalar first-order hyperbolic PDE. The recipe for writing the characteristic system can be found on pages 97-100 of [START_REF] Evans | Partial differential equations[END_REF]. In particular, the formula of interest to us is Formula (8) on page 98 of [START_REF] Evans | Partial differential equations[END_REF]. Specifically, since our equation is linear, the characteristic system becomes simpler (see Formula (17) on p.100 of [START_REF] Evans | Partial differential equations[END_REF]).

For the Vlasov equation under study, the following definition holds.

Definition 1. The set of characteristic curves of the linear gravitational Vlasov equation (1)-( 2) is the general solution of the following system of ordinary differential equations:

                 Ṁ (t) = 0, Ṫ (t) = 1, Ẋ(t) = V (t), V (t) = - ∇ x Φ(T (t), X(t)) M (t) . (7)
We introduce now the function ζ : R + → R, such that ζ(t) gives the solution f along the characteristic curves (M, T, X, V ) (i.e. ζ(t) = f (M (t), T (t), X(t), V (t))). On the characteristic curves, the Vlasov equation simply reduces to ζ(t) = 0.

The characteristic system takes hence the form

                       ζ(t) = 0, Ṁ (t) = 0, Ṫ (t) = 1, Ẋ(t) = V (t), V (t) = - ∇ x Φ(T (t), X(t)) M (t) , (8)
for t ∈ R + .

Proposition 1. For all s ∈ R + and for every (m, x, v) ∈ R * + × Ω t × R d , there exists a unique solution (M, T, X, V ) ∈ C ∞ of the initial value problem [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture[END_REF] with initial conditions

           M (s) = m, T (s) = s, X(s) = x, V (s) = v. (9)
Proof. It follows directly from the fact that the system ( 7) for (M, T, X, V ) satisfies the hypotheses of the Picard-Lindelöf Theorem and the sub-linear condition holds, since the points of singularity of Φ are excluded from the spatial domain Ω t for every time t ∈ R + .

Existence and uniqueness of the solution

We are now ready to employ what we have previously presented to obtain an existence and uniqueness result for System (1)-( 5)- [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF].

For the sake of simplicity, let us set z := (x, v), the initial data for the characteristic curves (X, V ); and Z(t; s, z) := (X, V )(t; s, x, v) the characteristic (X, V ) which start from (x, v) in the phase-space, i.e. it is equal to (x, v) when t = s, and which is parametrized by t.

The main tool to deal with the mobile domain in studying Equation ( 1)-( 2) is the backward absorbing time: Definition 2. The backward absorbing time τ Ω t (x, v) for a particle starting from x ∈ Ω t in the direction v ∈ R d , is defined as

τ Ω t (x, v) = inf{θ > 0 : X(θ; t, x, v) ∈ Γ t-θ }. If the set Θ := {θ > 0 : X(θ; t, x, v) ∈ Γ t-θ } is empty, then τ Ω t (x, v) = +∞.
In other words, τ Ω t (x, v) corresponds to the time of arrival at the border when we follow the characteristic X backwards from x ∈ Ω t with velocity v ∈ R d . This consideration follows from the two results summarized below: Proposition 2. The characteristic Z := (X, V ), solving (7), satisfies:

1. ∀t 1 , t 2 , t 3 ∈ R + and ∀z ∈ Ω t × R d , Z(t 3 ; t 2 , Z(t 2 ; t 1 , z)) = Z(t 3 ; t 1 , z); 2. ∀t, s ∈ R + , the map z → Z(t; s, z) is a C 1 -diffeomorphism and y → Z(s; t, y) (10)
is its inverse.

The proof of this result is classical and it will be omitted.

The infimum of θ > 0 such that X(θ; t, x, v) ∈ Γ t-θ is the time whereby, starting from it in y ∈ Ω θ , with velocity w ∈ R d , we arrive at the time t in X(t; θ, y, w) = x with a velocity v. So, if we are in x at the time t and we want to proceed backwards until we arrive in y at the time θ, we have to use the inverse function defines in (10):

x = X(t; θ, y, w) → X(θ; t, x, v), and it is precisely for this reason that, in Definition 2, we check X(θ; t, x, v).

The spatial domain Ω t of the Vlasov equation ( 1)-( 2) is not convex, for every t, so classical solutions of the boundary value problem for the Vlasov equation may not exist. Indeed, let us consider the characteristic part of the boundary, i.e.

Σ t 0 := N i=1 (x, v) ∈ ∂S t i × R d : (v -v i (x)) • n x = 0 (x, v) ∈ ∂P t × R d : (v -v P ) • n x = 0 .
Since Ω t is not convex, some velocity trajectories v from Σ t 0 can enter Ω t . It results that the method of characteristics -and therefore the explicit formula given born in the theorem above -does not define f on the points of Ω t interested by these trajectories.

However, Σ t 0 satisfies the hypotheses of Proposition 2.3 in [START_REF] Bardos | Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d'approximation; application à l'équation de transport[END_REF], and hence has zero Lebesgue measure.

We can express the solution of the initial-boundary value problem ( 1)-( 2)-( 5)-( 6) using the characteristics and the backward absorbing time: 2)-( 5)-( 6) associated to the initial condition f (m, 0, z) = f in (m, z), where z := (x, v). It is given by the same formula of the classical case, namely that, for a.e.

Theorem 1. If f in ∈ L 1 (R * + × Ω 0 × R d ) ∩ L ∞ (R * + × Ω 0 × R d ), then there exists a unique generalised solution f ∈ L 1 (R * + × R + × Ω t × R d ) ∩ L ∞ (R * + × R + × Ω t × R d ) of the boundary problem (1)-(
(m, t, z) ∈ R * + × R + × Ω t × R d , we have f (m, t, z) = f in (m, Z(0; t, z))1 {τ Ω t (z)>t} , (11)
where Z := (X, V ) represent the characteristic curve which solves [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture[END_REF].

Proof. The proof can be obtained by adapting, to external domains, the proof given by Bardos in [START_REF] Bardos | Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d'approximation; application à l'équation de transport[END_REF], which is based on semigroup theory. However, we give here a more direct proof, which is closer to the numerical strategy described in the next section, based on the study of the evolution of the unknown on the characteristic curves of the system. Let

N 0 ⊂ R * + × Ω 0 × R d be the set of zero Lebesgue measure in R * + ×R d ×R d such that f in is defined and of class L 1 ∩L ∞ on (R * + ×Ω 0 ×R d )\N 0 , and let N b ⊂ R * + × R + × Σ t -be the set of zero Lebesgue measure in R * + × R d × R d such that f is zero on (R * + × R + × Σ t -) \ N b . Consequently, the set N d := N b ∪ N 0 has zero Lebesgue measure in R * + × R d × R d .
We easily see that f , given by ( 11), satisfies the Vlasov equation almost everywhere. Indeed, Equation ( 11) defines a function on (R *

+ × R + × Ωt × R d ) \ N d such that f (m, t + s, Z(s; 0, z)) = f in (m, Z(0; t + s, Z(s; 0, z))) 1 {τ Ω t+s Z(s;0,z)>t+s} = f in (m, Z(0; t, z)) 1 {τ Ω t (z)>t} , for all s ∈ R such that Z(s; 0, z) ∈ Ω t × R d .
In fact, thanks to Proposition 2,

f in (m, Z(0; t + s, Z(s; 0, z))) = f in (m, Z(0; t, z))
because Z(0; t + s, Z(s; 0, z)) = Z(0; 2t, Z(2t; t + s, Z(s; 0, z))

Z(t;0,z):=y ) = Z(0; t, Z(t, 2t, y)) = Z(0; t, Z(t, 2t, Z(t; 0, z))) z .
Moreover,

1 {τ Ω t+s Z(s;0,z)>t+s} = 1 {τ Ω t (z)+s>t+s} = 1 {τ Ω t (z)>t} .
So, the function s → f (m, t + s, Z(s; 0, z)) is C 1 in the variable s for all (m, t, z) 11), we see moreover that

∈ (R * + × R + × Ωt × R d ) \ N d . Furthermore, we have that df ds (m, t + s, Z(s; 0, z)) = 0 for all (m, t, z) ∈ (R * + × R + × Ωt × R d ) \ N d and for all s such that (m, t + s, Z(s; 0, z)) ∈ R * + × R + × Ω t × R d . From Equation (
lim t→0 + f (m, t, Z(t; 0, z)) = lim t→0 + f in (m, Z(0; t, Z(t; 0, z))) = f in (m, z) for all (m, z) ∈ (R * + × Ω 0 × R d ) \ N 0 , whereas lim s→0 + f (m, t + s, Z(s; 0, z)) = f (m, t, z) = 0 for all (m, t, z) ∈ (R * + × R + × Σ t -) \ N b .
Therefore, f solves a.e. the Vlasov equation.

Uniqueness follows easily by noticing that if f is a generalised solution of Vlasov's equation, then the function

s → f (m, t + s, Z(s; 0, z)) is C 1 . Moreover, there exists N f ⊂ R * + × R + × Ω t × R d of zero Lebesgue measure such that df ds (m, t + s, Z(s; 0, z)) = 0 for all (m, t + s, Z(s; 0, z)) ∈ (R * + × R + × Ω t × R d ) \ N f
and for all s ∈] -min (t, τ Ω t (z)), 0[, and we have also

lim t→0 + f (m, t, Z(t; 0, z)) = f in (m, z) (m, z) ∈ (R * + × Ω 0 × R d ) \ N 0 ,
and lim

s→0 + f (m, t + s, Z(s; 0, z)) = 0 (m, t, z) ∈ (R * + × R + × Σ t -) \ N b .
Then, for all (m, t, z)

∈ (R * + × R + × Ω t × R d ) \ (N f ∪ N d ),
f (m, t, z) = lim ϵ→0 + f (m, 0, Z(0; t -ϵ, z)) = f in (m, Z(0; t, z))
and for t > τ Ω t (z), f (m, t, z) = 0.

Description of the numerical strategy

In this section we describe the numerical method used for the simulation of equations ( 1)-( 5)-( 6) in two dimensions in space and velocity. Our procedure, based on a particle method and a splitting strategy, is an evolution of the approach introduced in [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF]. The substantial difference of our problem from the one studied in [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF] is that the dust particles in planetary rings, unlike gas molecules, being affected by the gravitational acceleration due to the planet and moons, satisfy the Vlasov equation ( 1), which is more complicated to deal with than the free transport equation. On the other hand, the boundary condition are, in this problem, particularly easy to handle, because the surfaces of the satellites are modelled with an absorbing boundary: the particles that collide with it disappear from the problem's domain. We place ourselves in the reference frame of the planet, which is thus supposed at rest. The centres of mass of the moons evolve independently according to the following equations in polar coordinates:

ξ i (t) = r i (t) cos θ i (t) sin θ i (t) with r i (t) = r 0 i θ i (t) = ω i t , for i = 1, . . . , N . (12)
In [START_REF] Evans | Partial differential equations[END_REF], ω i is the constant angular velocity of the circular orbit of the i-th moon and r 0 i is its constant distance from the origin of the coordinate system. If the i-th moon has an elliptical shape, then we have to consider also its rotation around its own centre of mass. So, if δ i describes the angle formed by the major axis of the i-th moon and the x-axis of the coordinate system, and if σ i represents the angular velocity of the i-th satellite around its own centre, we have δ i (t) = δ i (0) + σ i t.

The initial density of the dust particles is discretized by means of a collection of weighted smooth shape functions centred on the numerical particle positions and velocities (x 0 k , v 0 k ) 1⩽k⩽N 0 , that is

f in N 0 ,ε (x, v) = N 0 k=1 α k φ ε (x -x 0 k ) φ ε (v -v 0 k ), (13)
where N 0 represents the initial number of numerical particles, α k is the weight of the k-th numerical particle (which represents α k "real" particles). In [START_REF] Garber | Subjects great and small: Maxwell on Saturn's rings and kinetic theory[END_REF], the shape function φ ε (x) = φ(ε -1 x)/ε d is a smooth function with compact support. The term "numerical particles" is here used for avoiding any confusion with the (real) number of dust particles. Once the number N 0 of numerical particles has been chosen, the initial positions (x 0 k ) 1⩽k⩽N 0 and velocities (v 0 k ) 1⩽k⩽N 0 are sampled according to the initial density f in thanks to a Monte-Carlo procedure. Then, the positions and velocities of the numerical particles evolve in time, according to the explicit Euler scheme, with time step ∆t

   x n+1 k = x n k + ∆tv n k v n+1 k = v n k -∆t ∇xΦ(t n ,x n k ) m 1 ⩽ k ⩽ N n , (14)
where N n is the number of numerical particles at the time t n (which can differ from N 0 because some particles can be absorbed by the planet or the moons between t 0 and t n ).

In order to simplify the interactions between the moons and the numerical particles, we use a time-splitting between the transport (free flow of the particles in the absence of any interaction, mathematically represented by the transport operator v •∇) and the treatment of the differents interactions (flow of the particles due to the planet-moons gravitational field, and the absorbing boundary condition on Σ t -). In other words, in each time interval [t n , t n+1 [, we first move the moons and the planet independently of the motion of the dust particles, then we freeze them and we transport the numerical particles and perform the treatment of the boundary conditions. Thus, when we move the numerical particles, the domain is fixed, allowing us to come back to deal with the boundary conditions of a fixed domain instead of working in a mobile domain. Once determined the positions of the largest bodies at the time t n+1 , we move the N n macro-particles, according to the equations of motion [START_REF] Goldreich | The dynamics of planetary rings[END_REF], finding (x n+1 k , v n+1 k ) 1⩽k⩽N n . Then, we test on every numerical particle k if x n+1 k ∈ Ω n+1 and otherwise we impose the boundary condition. To do so, we only need a cartesian equation of the surface of the planet and of the satellites, which for elliptical bodies with axis lengths a and b is

(x, y) ∈ R 2 | (x -r i (t) cos(θ i (t)) a 2 + (y -r i (t) sin(θ i (t)) b 2 ⩽ 1 .
Once written the coordinates of macro particles in the reference frame of the ellipse, we control if they verify the cartesian equation or not. If yes, we apply the boundary condition, namely we remove the k-th particle from the domain, then we renumber the remaining numerical particles. The number of numerical particles at time t n+1 hence becomes then N n+1 ⩽ N n . We neglect the mass increase of the moon coming from the absorption of particles.

The shortcoming of this strategy is that, when a particle collides with the i-th body, it does not allow us to determine the position of the intersection between the particle's trajectory and the i-th moon. However, since we are dealing with absorption boundary conditions, the value of this intersection becomes superfluous for our discussion. Although this strategy is less precise than the one in which the intersection is calculated, the graphical results of these two approaches are similar, with the advantage that the former allows us to deal with bodies having complicated shapes and is less computationally expensive.

Another problem that could arise is shown in the Figure 1: during the time interval [t n , t n+1 [ a dust particle collides with a moon, but the method does not detect the collision because the particle is outside the moon at time t n+1 . However, this problem can be controlled by reducing the time step ∆t. The value of ∆t should be chosen globally for all numerical particles in order to guarantee i) the partial superposition between the domain representing the moon at time t n and at time t n+1 and ii) if R > 0 is the radius of the smallest moon, max k=1,...,

N n ∥x n+1 k -x n k ∥ < 2R.
If the previous constraint is not satisfied for some t n , the time step ∆t has to be reduced (for all particles) in order to fullfill it. The computational cost of our method is clearly proportional to N 0 /∆t, because there are no interactions between the numerical particles.

Figure 1: A situation in which the k-th particle collides with the i-th moon, but the numerical method does not notice the collision.

Numerical results

We now describe some simulations implemented for the planetary rings problem in two spatial dimensions (i.e. four dimensions in phase-space). The data used in the numerical simulations are scaled in order to be compatible with the data of the giant planets of the Solar system. Here we have chosen moons of sufficiently high size and mass in order to have appreciable and visible results.

Simulation 1

In this scenario, we simulate the orbit of an elliptical shepherd moon around a planet and observe how the moon's motion creates a separation in the dust cloud around the planet, thus contributing to the formation of a ring system. We consider a circular planet, with radius ρ P and mass m P , and an elliptical shepherd moon, with axes a 1 and b 1 and mass m S 1 . We moreover assume that the moon, with center of mass ξ 0 = r 0 1 cos θ 0 1 sin θ 0 1 at the time t = 0, orbits the planet with constant angular velocity ω 1 and rotates around its own centre of mass with constant angular velocity σ 1 . The quantity δ 0 1 is the initial angle from the semi-major axis of the i-moon and the x-axis of the coordinate system.

The spatial domain of the planetary ring problem is, for d = 2:

Ω t := R 2 \ P t ∪ S t 1 , t ⩾ 0.
Our numerical spatial domain is a truncation of Ω t 1 for each time t n :

Ωn := D \ P n ∪ S n 1 ( 15
)
where We have supposed that the initial distribution f in is factorized as the product of a function depending on the space variable only and a function depending on a suitably chosen velocity v orb = v orb (x). More precisely, let r min = 6.7 × 10 7 m, r max = 1.8 × 10 8 m and define

D := [-1, 1] × [-1,
R := (x 1 , x 2 ) ∈ D, r min ⩽ x 2 1 + x 2 2 ⩽ r max and S 2 orb (x) := v ∈ R 2 , ∥v∥ = v orb (x) . The initial distribution is f in (x, v) = 1 R (x)1 S 2 orb (x) (v). ( 16
)
We thus first consider the finite set (x 0 k ) 1⩽k≤N 0 as a realisation of the density function associated to 1 R (x): we generate the positions (x 0 k ) 1⩽k⩽N 0 as the realisations of a uniform density on R, in a probabilistic way. The initial velocity v 0 k of the k-th macro-particle is then chosen as the orbital velocity of a point at x 0 k (i.e. we suppose that the radial velocity v r is zero). Thus, written its initial position

x 0 k = r 0 k cos Θ 0 k sin Θ 0 k in polar coordinates, its
orbital velocity is:

v 0 k := v k,orb (-sin Θ 0 k , cos Θ 0 k ), k = 1, . . . , N 0 , (17)
where

v k,orb = Gm P r 0 k (18)
is the approximation of the modulus of the orbital velocity for dust particles. The particle weights are identical, i.e. we suppose that

α k := ||f in || L 1 (Ω 0 ×R 2 ) N 0 = n 0 |R| N 0 , 1 ⩽ k ⩽ N 0 , (19)
where n 0 :=

||f in || L 1 (Ω 0 ×R 2 ) |R|
. This quantity has been normalized to 1 in the simulations.

Then, for every n, we move the corresponding numerical particle in the time interval [t n , t n+1 ] using an explicit Euler method and we check if the numerical particle falls into the planet or the moons. If so, it is eliminated.

The next step consists in the reconstruction of the density f n+1 N 0 ,ε . As in [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF], we use Bsplines of 3-order in two space dimensions as shape functions ϕ ε , with shape sizes

ε 1 = h 1/2 1 and ε 2 = h 1/2
2 , where h 1 and h 2 are the initial distances between two numerical particles in the first and second direction, respectively.

In order to visualize the space-time evolution of the dusty particles, we reconstruct the spatial density

ρ(t, x) := R 2 f (t, x, v)dv
from the positions and the velocities of the numerical particles and we plot the numerical macroscopic density

ρ n N 0 ,ε (x) := R 2 f n N 0 ,ε (x, v)dv = N n k=1 α k ϕ ε (x -x n k ).
In Figure 2 we show the results of the simulation. We have used N 0 = 587 116 initial numerical particles. In Saturn rings, some estimates suggest that the order of magnitude of the number of dust particles is 3 × 10 16 . According to [START_REF] Kawai | An inelastic collision model for the evolution of "planetary rings[END_REF], it means that a numerical particle corresponds to 5 × 10 10 dust particles. The planet and the moons are coloured in red to facilitate their visualization. One can notice that the moons sweep away the dust in front of them as they orbit the planet, creating a path where the density ρ n N 0 ,ε of the rings is very low. This result is consistent with astronomers' assumptions about the creation of ringlets based on data collected by the Cassini spacecraft [START_REF] Horányi | Diffuse Rings[END_REF].

Simulation 2

In this section we describe the results of the simulation of a close-up view of the moon in a previously open gap, neglecting the planet's gravitational role on the moon and ring dust. Let us then consider the square D and take as domain, at each discrete time t n , the set: Ωn := D \ S n . The space in which the dust particles are initially located is the following:

R := {(x 1 , x 2 ) ∈ D | r min ⩽ x 2 ⩽ r max }. ( 20
)
where r min = -5 × 10 7 m and r max = 5 × 10 7 m. As initial distribution, we consider which is uniform in space and non-zero only if the velocity is zero in each component. So, (x 0 k ) 1⩽N 0 are selected as realisation of a uniform density, while (v 0 k ) 1⩽N 0 are set to zero. The weights (α k ) 1⩽k⩽N 0 are chosen as in [START_REF] Kawai | An inelastic collision model for the evolution of "planetary rings[END_REF], with R defined in [START_REF] Kenworthy | Mass and period limits on the ringed companion transiting the young star J1407[END_REF].

f in (x, v) = 1 {x∈R} δ 0 (v), (21)
The physical data for the moon are the following:

Moon m S = 6.0 × 10 18 kg a = 1.50 × 10 6 m b = 1.00 × 10 6 m δ 0 = 3.14/3 rad

The moon is supposed to be immobile, in order to study only the influence of its gravitational force on the particles' ring. At each time step, we move the particles, following the numerical solution of the equation of motion [START_REF] Goldreich | The dynamics of planetary rings[END_REF] in the time interval ]t n , t n+1 ]. Particles which fall into the moon are eliminated. Then, we reconstruct the density f n+1 N 0 ,ε as described in Section 5.1 and we employ it to define the spatial numerical density

ρ n+1 N 0 ,ε (x) = N n+1 k=1 α k ϕ ε (x -x n+1 k ).
Figure 3 shows the time evolution of the numerical density ρ n N 0 ,ε on D. In this simulation, we have employed 485 825 numerical particles.

By comparing the two pictures, we observe the effect of the gravitational attraction by the moon on the dust particles and, moreover, the emergence of wavy, jagged edges, consistent with what has been observed in Saturn's rings [START_REF] Esposito | Planetary rings: a post-equinox view[END_REF].

This effect is not due to the particle method, because the stochasticity appears only in the discretization of the initial condition. The other steps of our strategy are deterministic and therefore the oscillations are a consequence of the model itself.

Conclusion

With this paper, we aimed to provide a robust mathematical model capable of describing the formation of gaps and patterns within planetary ring systems. Firstly, we focused on a formal justification of the proposed Vlasov model ( 5)- [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF]. The physical processes involved in a ring system are numerous, yet some can, at first glance, be neglected. A careful selection of those relevant to the problem of our interest was therefore necessary. Subsequently, we proved, via the method of characteristics on a time-dependent domain, the well-posedness of the Vlasov equation ( 5)- [START_REF] Charles | Mathematical and numerical study of a dusty Knudsen gas mixture: Extension to non-spherical dust particles[END_REF]. The proof of Theorem 1 via this method offers a natural connection to the presented numerical strategy. Indeed, ultimately, the purpose of the paper was to provide a numerical method to simulate the model. We used the particle method coupled with a splitting technique, useful for capturing collisions between dust particles with the planet or moons. The simulations that emerge (Figures 2 and3) are coherent with the predicted theoretical results.

In future works, collisions between particles and certain source terms (such as the presence of a moon, e.g. Titan, capable of supplying the rings with new particles) can be included and thus more advanced models can be studied.

m r 0 1 = 1 . 2 × 10 8 m θ 0 1

 1121 1] (scaled according to the length L 0 = 2 × 10 8 m) is the square domain for the simulations.The table shows the values of the data we have selected for these simulations.PlanetMoon m P = 5.70 × 10 26 kg m S 1 = 8 × 10 21 kg ρ P = 5.82 × 10 7 m a 1 = 1.6 × 10 6 m b 1 = 1.6 × 10 6 = 0.00 rad ω 1 = 1.33 × 10 -3 rad/s δ 0 1 = 1.05 rad σ 1 = 4 × 10 -10 rad/s

Figure 2 :

 2 Figure 2: Time history of the spatial numerical density ρ n N 0 ,ε at instants t = 0 s, t = 6.0 × 10 5 s, t = 1.2 × 10 6 s, t = 1.8 × 10 6 s, starting from the image at the top left, proceeding from left to right and from top to bottom. The axis are scaled according to the length L 0 = 2 × 10 8 m.

Figure 3 :

 3 Figure 3: Time history of the spatial numerical density ρ n N n ,ε at instants t = 0 s (left) and t = 4.96 × 10 4 s (right). The axis are scaled according to the length L 0 = 2 × 10 8 m.

Table 1 :

 1 Data of Solar system planets with planetary rings.

		Saturn	Uranus	Neptune	Jupiter
	m P	5.7 × 10 26 kg	8.7 × 10 25 kg	1.0 × 10 26 kg	1.9 × 10 27 kg
	d P	1.2 × 10 8 m	3.8 × 10 7 m	6.3 × 10 7 m	1.2 × 10 8 m
	d S	1.4 × 10 12 m	2.7 × 10 12 m	4.5 × 10 12 m	7.8 × 10 11 m

Acknowledgments: Work supported by the Austrian Exchange Service ( ÖAD) in the framework of the Austrian-French projects, grants FR 01/2021 and MULT 11/2020. The second author acknowledges partial support from the Austrian Science Fund (FWF), grants P33010, W1245, and F65. Article written under the auspices of the Italian National Institute of Higher Mathematics (INdAM), GNFM group.

Conflicts of Interest Statement:

The authors confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.