Silvia M Lenzi

Frédéric Nowacki

Alfredo Poves

Kamila Sieja

Kamila Sieja The

The land of deformation south of 68 Ni

Keywords: numbers: 21.60.Cs, 21.10.-k, 21.10.Re

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

I. INTRODUCTION

In the last decades, more and more experimental evidence has been accumulated establishing the breaking of the shell closures known at the stability valley when approaching the drip lines, mainly at the very neutron-rich side. The first example of an unexpected disappearance of a shell closure (N = 8) was found in 11 Be, whose ground state is an intruder 1/2 + located 320 keV below the "natural" 0 ω 1/2 -state [1]. However, the true relevance of this finding was not recognized till many years later, and even now it is shadowed by the fame of its more neutron-rich isobar 11 Li. In the sd shell, the expected N = 20 semi-magic neutron rich nuclei turned out to be actually well deformed [2,3]. The presence of deformed ground states in the nuclei around 32 Mg and the breaking of the N = 20 shell closure in this region has been extensively studied experimentally and theoretically (see [4][5][6] for shell model reviews). The aim of these studies is to map the limits of the so-called island of inversion, i.e. the region of nuclei where the strong quadrupole correlations overcome the spherical mean-field gaps, favoring energetically the deformed intruders, which often become ground states. Indeed, another basic aim is to understand microscopically the dynamics responsible for these shape transitions.

The question of the persistence of the N = 40 harmonic oscillator closure in neutron-rich nuclei comes in naturally in this context. The high lying 2 + state observed in 68 Ni and its low B(E2; 2 + → 0 + ) value are the result of the relatively large energy gap separating the pf and 0g 9/2 orbitals [7]. However, this gap gets reduced (or even disappears) when protons are removed from 68 Ni: The nucleus 66 Fe, with only two protons less, shows a sudden change in nuclear structure with an increased collectivity manifested via its very low lying 2 + state. Along the iron chain, indications for a collective behavior come from the systematics of the 2 + states [8] as well as from the recent measurement of the B(E2) values in 64,66 Fe [9,10]. The evolution of the B(E2) values in iron isotopes points to a sudden increase of collectivity when approaching N = 40. Only very recently the first measurement of the excited levels in 64 Cr has been reported [11]. This is the nucleus where theoretical calculations predicted the lowest lying 2 + level in the region [12][13][14]. The measured 2 + state energy agrees within 100 keV with those theoretical predictions. The sudden shape change at N = 40 challenges the theoretical models, which have to account for a particularly rapid shell evolution responsible for these effects. Recent beyond mean-field HFB+GCM calculations with the Gogny force, reported in Ref. [15], show an increase of collectivity towards the proton drip line. The spherical neutron single-particle energies obtained in the HFB approach reveal almost no variation of the N = 40 gap with the proton number between Z = 20 and Z = 28. As a consequence, only a moderate collectivity in the iron and chromium chains is found, without visible structure changes between them. The collectivity at N = 40 has been previously a subject of many shell model studies using different valence spaces and interactions [11][12][13][14]16]. In particular, it has been shown that the shell model calculations using the 0f 5/2 , 1p 3/2 , 1p 1/2 , 0g 9/2 neutron orbits (the f pg valence space) and realistic interactions [13], can reproduce rather well the level schemes of 62,64 Fe, but fail to do so for the 2 + state of 66 Fe [16]. To reproduce the large quadrupole collectivity in this mass region, the inclusion of the neutron 1d 5/2 orbital is needed, as first surmised in Ref. [12] and confirmed recently in Ref. [9]. This can be explained in terms of the quasi-SU3 approximate symmetry: In this framework, the deformation is generated by the interplay between the quadrupole force and the central field in the subspace consisting on the lowest ∆j = 2 orbitals of a major shell [17].

In this work we discuss in detail how the sudden onset of collectivity can be interpreted in terms of shell model calculations in large model spaces and we look for the similarities in the deformation-driving mechanisms at N = 20 and N = 40 when approaching the neutron drip line. We present novel calculations for the N = 40 region, in a model space including the pf -shell for protons and the 0f 5/2 , 1p 3/2 , 1p 1/2 , 0g 9/2 and 1d 5/2 orbits for neutrons. We discuss first our model space and interaction as well as the computational details in Sec. II. The shape evolution along the N = 40 line is studied in Section III. Next, we illustrate the development of deformation along the isotopic chains of iron and chromium in Sec. IV and we discuss the structure of the nickel isotopes. Conclusions are given in Sec. V.

II. THE VALENCE SPACE

The physics of the N = 40 nuclei has been discussed in the recent past both experimentally and theoretically [5,8,9,[11][12][13][14][15][16]. In 68 Ni, the N = 40 harmonic oscillator shell closure is somewhat weakened, however, whereas some of its properties seem consistent with a superfluid behavior, others may point to a double magic character. Below 68 Ni, the iron and chromium isotopic chains are characterized by an open proton shell which favors the development of quadrupole correlations. Such a situation occurs also for the N = 20 nuclei in the so called island of inversion around 32 Mg [18]: the protons occupy the N = 2 harmonic oscillator shell and the neutrons may lie either in the N = 2 harmonic oscillator shell (normal filling) or in the 1p 3/2 and 0f 7/2 orbitals from the N = 3 major shell (intruder configurations). The latter orbits form a quasi-SU3 block [19] which enhances the quadrupole correlations. Extended shell model calculations by the Tokyo group, with the SDPF-M effective interaction [18], predict the dominance of different type of configurations (0p0h, 1p1h, 2p2h) for the ground states of the nuclei with 18 ≤ N ≤ 22 and 10 ≤ Z ≤ 14, delineating the contours of the island of inversion.

To assess if such a scenario could develop around 68 Ni, we adopt here a model space based on a 48 Ca core which comprises the pf -shell for protons and the 1p 3/2 , 1p 1/2 , 0f 5/2 , 0g 9/2 , and 1d 5/2 orbits for neutrons. The degrees of freedom that can be encompassed in this valence space are very similar to the ones present in the study recalled previously [18], and the same kind of phenomena could therefore be described. The advantage of this valence space is that it contains all the physical degrees of freedom important for the description of the low-lying properties of these nuclei, it is computationally tractable, (see precisions hereafter) and it is almost free of center of mass spuriousity since its main components, the 0f 7/2 → 0g 9/2 excitations, are excluded from the space.

The interaction proposed in this work, denoted hereafter LNPS, is a hybrid one, based on several sets of realistic two-body matrix elements (TBME). Its main building blocks are:

• The last evolution of the Kuo-Brown interaction (KB3gr) for the pf -shell [20];

• The renormalized G-matrix of Ref. [21] with the monopole corrections introduced in [22], for the remaining matrix elements involving the 1p 3/2 , 1p 1/2 , 0f 5/2 , and 0g 9/2 neutron orbits.

• The G-matrix based on the Kahana-Lee-Scott potential [23], for the matrix elements involving the 1d 5/2 orbit. This potential has been successfully employed in the definition of the recent SDPF-U shell model interaction [24] for the description of neutron-rich sdpf nuclei.

In addition, another set of experimental constraints has been taken into account for the final tuning of the monopole Hamiltonian:

• The Z = 28 proton gap around 78 Ni is inferred from recent experimental data in 80 Zn and fixed so as to reproduce the measured B(E2). A standard polarization charge of 0.5e has been used in all the calculations presented in this work.

• The size of the N = 50 neutron gap in 78 Ni has been estimated to be ∼ 5 MeV. The evolution of the neutron gap with the neutron number is rather independent on the proton number. On the contrary, the systematics [25] shows that the 0g 9/2 -1d 5/2 gap in Zr isotopes increases by 3 MeV when the neutron 0g 9/2 orbital is filled. We assume a similar behavior for the Ni chain. With this assumption, the observed 5 + , 6 + states in 82 Ge, which are supposed to be 1p-1h excitations across the N = 50 gap, are correctly reproduced [26].

Finally, in order to compensate for the absence in our space of the third member of the quasi-SU3 sequence, the 2s 1/2 neutron orbit, we have increased the quadrupolequadrupole interaction of the neutron 0g 9/2 and 1d 5/2 orbits by 20%. The effective single particle energies (ESPE) of the LNPS interaction are shown in Fig. 1 for the neutron orbits at N = 40, between 60 Ca and 72 Ge. For comparison, we show in the same figure the neutron ESPE at N = 20 between 28 O and 36 S, calculated with the SDPF-U interaction [24]. The similarities are striking. In the N = 20 case, a reduction of the neutron 0d 3/2 -0f 7/2 gap takes place when protons are removed from the proton 0d 5/2 orbital. This feature, accompanied by the proximity of the quadrupole partner neutron orbitals 0f is responsible for the formation of the island of inversion at N = 20. At N = 40 one observes the same behavior but for the neutron 0f 5/2 -0g 9/2 gap when going down from 68 Ni and the closeness of the quadrupole partners 0g 9/2 -1d 5/2 .

The analogy in the shell evolution between N = 20 and N = 40 suggests therefore a possibility of another island of inversion below 68 Ni.

As mentioned already, the deformation driving role of the neutron 1d 5/2 orbital below 68 Ni was already discussed in the shell model calculations of Ref. [12]. However, due to the computational limitations at that time, the calculations were performed with a closed neutron 1p 3/2 orbit, therefore using 52 Ca as a core. This is no longer an issue here: The present calculations have been carried out including up to 14 particle-14 hole excitations across the Z = 28 and N = 40 gaps, when it appeared necessary to assure the convergence of the calculated electromagnetic properties. The largest dimensions of the matrices treated here reach 10 10 in the case of 64 Fe. All the calculations of this work have been performed using the m-scheme shell model code ANTOINE [27]. The change of structure in even-even N = 40 isotones from 68 Ni down to 60 Ca is illustrated in Fig. 2. In part (a) the evolution of the excitation energy of the 2 + state is plotted while part (b) shows the corresponding evolution of the E2 transition rates. The available experimental data are satisfactorily reproduced in all cases. Our theoretical approach predicts the maximal deformation at the middle of the proton shell, i.e. in 64 Cr, where the calculated 2 + energy is the lowest and the value of the B(E2; 2 + → 0 + ) is the largest.

In Table I we list the extra occupancies of the two neutron intruder orbitals, 0g 9/2 and 1d 5/2 , relative to the normal filling. The occupation of the 0g 9/2 orbit in 68 Ni is close to 1, which corresponds to ≈50% of the calculated wave function having a J = 0 pair excited across N = 40 gap. The wave function is however difficult to interpret, because, in spite of its very large content of excited pairs, the doubly magic component is still substantial. This is probably the reason why 68 Ni shows at the same time an increase of the 2 + excitation energy typical for doubly magic nuclei, and no sign of shell closure in the neutron separation energy.

The occupation of both neutron intruder orbits grows rapidly when protons are removed, due to the reduction of the neutron N = 40 gap shown in Fig. 1. The ratio of the 0g 9/2 to 1d 5/2 occupations also evolves from Ni to Ca, partly as an effect of the level crossing taking place around 62 Ti, but mostly due to the quasi-SU3 structure of the intruder states. The very important role of the 1d 5/2 neutron orbit in the build up of collectivity in this region sheds doubts about the results of Ref. [14], which do not include this orbit. The more so when one examines the very unrealistic ESPE that they enforce into their schematic interaction in order to get results close to the experimental data.

In Table I we list as well the percentages of neutron n-particle n-hole excitations in the ground state wave functions. The 4p-4h components are dominant in Fe, Cr, Ti and Ca, however 2p-2h and 6p-6h contributions are sizable. The complexity of the wave functions constitutes the main difference between N = 40 and N = 20 regions. In the latter, nearly pure 2p-2h components has been shown to dominate the 0 + ground states [18].

It should be also pointed out that this evolution of the neutron filling and of the particle-hole structure does not mean that the nuclei will become more and more deformed with decreasing Z; the ground-state deformation properties result from the total balance between the monopole and the correlation energies (mainly of a proton-neutron character). In Table I, we list also these correlation energies extracted from the multipole Hamiltonian along the N = 40 line. Indeed, the correlation energy, reflecting the deformation, increases from Ni to Cr, where it reaches its maximum, and then diminishes toward Ca. The transition between Ni and Cr is not gradual: The removal of two protons already provokes and abrupt change, from spherical to a strongly deformed prolate shape.

IV. THE EVOLUTION OF THE DEFORMATION ALONG THE ISOTOPIC CHAINS

The isotope chains provide a very illustrative picture of the evolution of deformation in this region. In Fig. 3 we show the results obtained for the iron chain. In the panel (a) the excitation energies of the 2 + states are compared to the available experimental data. In panel (b) we show the ratio of the excitation energies E(4 + )/E(2 + ) and E(6 + )/E(4 + ), which make it possible to recognize whether a nuclear spectrum is close to that of the rigid rotor or not. Let us remind that the perfect rotational limit would require E(4 + )/E(2 + )=3.33 and E(6 + )/E(4 + )=2.1. Panel (c) depicts the theoretical and experimental B(E2) transition probabilities along the yrast bands. The measured B(E2; 2 + → 0 + ) values from Ref. [9] are plotted with dots and from Ref. [10] with crosses. Finally, the intrinsic quadrupole moments derived from the calculated spectroscopic ones are shown in part (d). To establish a connection between the laboratory and the intrinsic frames we use the relations:

Q int = (J + 1)(2J + 3) 3K 2 -J(J + 1) Q spec (J), K = 1, (1)
and

B(E2, J → J -2) = 5 16 e 2 | JK20|J -2, K | 2 Q 2 int (2)
for K = 1 2 , 1. The same properties are plotted in Fig. 4 for the chromium chain. An excellent agreement with the experiment is found for all excitation energies.

The known B(E2) transition rates are well reproduced within the error bars as well. Comparing the results for both isotopic chains it can be seen that the onset of deformation occurs at different neutron number. In the chromium chain the intrinsic quadrupole moment stays constant along the yrast the band already at N = 38. This is one of the fingerprints of a good rotor behavior. The iron isotopes undergo the transition at N = 40. We have also verified that the intrinsic quadrupole moments obtained from the spectroscopic ones (Eq. ( 1)) are nearly equal to those obtained from the transition probabilities according to Eq. ( 2) in 62,64,66 Cr and 66,68 Fe. We obtain a value of Q int ∼ 150efm 2 in 62,64,66 Cr which corresponds to β ∼ 0.35. The deformation of 66,68 Fe is slightly lower, with Q int = 145efm 2 (β ∼ 0.3). In the lighter Fe and Cr nuclei, the discrepancies between the calculated quadrupole moments and those obtained in the rotational scheme from the B(E2)'s are much larger.

Let us further discuss the structure of the calculated states. In Table II we show the occupancies of the neutron intruder orbitals in the ground states of the chromium and iron chains. In the case of the 0g 9/2 orbit the extra occupancy is reported, i.e. the difference between the value obtained in the configuration mixing calculation and the one that corresponds to the normal filling. The neutron intruder occupations increase with N in both isotopic chains, however the absolute occupancies are larger in the chromium chain. As mentioned already for the N = 40 nuclei, the strong deformation in the Cr chain is not only due to the increased population of the 0g 9/2 -1d 5/2 doublet, but to the strong protonneutron correlations which tend to be maximal when four protons are active in the pf shell as well.

In Fig. 5 we show also the evolution of the collectivity in the nickel chain: in (a) the theoretical energies of the first excited 2 + states are compared with experimental ones, while in (b) we show the calculated transition rates in comparison to the available data. The agreement in the calculated energies is very good for all nickels. Concerning transition rates, the model reproduces well the systematics with the minimum in 68 Ni and the rapid increase of collectivity in 70 Ni. However, the calculated value is closer to the lower tip of the error bar. The transition rate in 64 Ni is underestimated. In this case we know that a better agreement with experiment can be obtained in a full pf shell calculation as the neutron excitations from 0f 7/2 orbital are here more important than those through the N = 50 gap and, as expected, the occupation of the 1d 5/2 orbit remains close to zero in the calculations.

Finally, let us note the recent work on 68 Ni [28], where the excited 0 + states have been investigated and a candidate for a proton 2p-2h intruder has been proposed at an excitation energy of 2.2 MeV. We have calculated the excited 0 + states in our model space and we obtain the first excited 0 + 2 state at an energy of 1.2 MeV and the second one 0 + 3 at an energy of 2.4 MeV, in good agreement with experiment. The too low excitation energy of the first state can be a result of the delicate mixing between the 0 + 1 ground state and the excited 0 + 2 state. These two states are characterized by similar proton occupancies with leading 0p-0h (neutron) configuration for the 0 + 1 ground state and 2p-2h (neutron) configurations for the 0 + 2 . For the 0 + 3 state, located at 2.4 MeV, the striking feature is that the dominant proton configuration has exactly two 0f 7/2 protons less than the ground state, i.e. a pure 2p-2h proton configuration.

The total quadrupole sum rules for these states amount to 904, 1162 and 2025 e 2 fm 4 for the ground state, the 0 + 2 , and the 0 + 3 states respectively. This shows in particular that the 0 + 3 state carries moderate deformation, despite of its proton intruder nature.

V. CONCLUSIONS

The aim of this work was to study the rapid onset of collectivity below 68 Ni, suggested by the experimental evidence, which is a great challenge for any theoretical model. Here we have presented shell model calculations in a large valence space including the pf shell for protons and the 0f 5/2 , 1p 3/2 , 1p 1/2 , 0g 9/2 , and 1d 5/2 orbitals for neutrons. The effective interaction for this model space has been build up using sets of realistic TBME as a starting point and applying monopole corrections thereafter. The progresses in algorithms and computer power have made it possible to achieve the largest shell model diagonalizations in this region of nuclei up to date.

The present calculations have been performed in Ni, Cr, Fe, Ti and Ca isotopes around N = 40. We have found a satisfactory agreement between the theoretical results and the available data, for both excitation energies and transition rates. In particular, it has been pos-sible to describe correctly the rapid onset of collectivity below 68 Ni without any collapse of the spherical gaps.

It has been also shown that the onset of deformation develops at N = 40 in the iron chain and already at N = 38 in chromium isotopes. The maximum deformation is predicted for the chromiums which exhibit features typical of rotational nuclei. The continuous advances in the experimental side, with particular regard to the transition probability measurements, will be a stringent test for these theoretical predictions.

The calculated wave functions of the deformed ground states were shown to contain large amounts of many particle-many hole configurations, with around four neu-trons occupying intruder orbitals 0g 9/2 , 1d 5/2 . We have also related the observed shape change with the underlying evolution of the spherical mean-field, which bears many similarities with the one at the N = 20 island of inversion.

2 FIG. 1 :

 21 FIG. 1: (Color online) Neutron effective single particle energies obtained with the LNPS interaction at N = 40 (upper part) and with the SDPF-U interaction at N = 20 (lower part).

  III. THE PROPERTIES OF THE N = 40 ISOTONES FOR Z≤28 AND THE NEW REGION OF DEFORMATION

FIG. 2 :

 2 FIG. 2: (Color online) Evolution of the nuclear structure along the N = 40 isotonic chain: (a) theoretical 2 + excitation energies and (b) B(E2; 2 + → 0 + ) values, compared to the available experimental data.

FIG. 3 :

 3 FIG. 3: (Color online) Theoretical results along the iron isotopic chain in comparison with the available experimental data: the excitation energies of the 2 + states are shown in panel (a), in panel (b) we present the ratio of energies of E(J +2)/E(J), the B(E2) transition rates are plotted in panel (c) and the calculated intrinsic quadrupole moments in panel (d). Two experimental sets of the B(E2) values are shown: from Ref. [9] (black dots) and from Ref. [10] (crosses).

FIG. 4 :

 4 FIG. 4: (Color online)The same as in Fig.3but for the chromium isotopes.

FIG. 5 :

 5 FIG. 5: (Color online) Theoretical versus experimental energies and transition rates of the nickel isotopes in the vicinity of N = 40.

TABLE I :

 I Occupation of the neutron intruder orbitals and percentage of particle-hole excitations across the N = 40 gap in the ground states of the N = 40 isotones. The last column contains the correlation energies evaluated for these states.

	Nucleus νg 9/2 νd 5/2 0p0h 2p2h 4p4h 6p6h Ecorr
	68 Ni 0.98 0.10 55.5 35.5 8.5 0.5 -9.03
	66 Fe	3.17 0.46	1	19	72	8 -23.96
	64 Cr 3.41 0.76	0	9	73	18 -24.83
	62 Ti	3.17 1.09	1	14	63	22 -19.62
	60 Ca 2.55 1.52	1	18	59	22 -12.09

TABLE II :

 II Extra occupations of neutron 0g 9/2 and 1d 5/2 orbitals in the ground states of the chromium and iron chains.

	Nucleus N ν0g 9/2 ν1d 5/2
	62 Fe 36 0.95	0.12
	64 Fe 38 2.0	0.27
	66 Fe 40 3.22	0.51
	68 Fe 42 2.30	0.62
	60 Cr 36 1.55	0.31
	62 Cr 38 2.77	0.66
	64 Cr 40 3.41	0.76
	66 Cr 42 2.28	0.90

Acknowledgments. This work is partly supported by a grant of the Spanish Ministry of Science and Innovation MICINN (FPA2009-13377), by the IN2P3(France)-CICyT(Spain) collaboration agreements, by the Spanish Consolider-Ingenio 2010 Program CPAN (CSD2007-00042) and by the Comunidad de Madrid (Spain), project HEPHACOS S2009/ESP-1473.


Philippe Robutel

Nicolas Rambaux

Julie Castillo-Rogez

Analytical description of physical librations of Saturnian coorbital satellites Janus and Epimetheus

) have proposed a model of rotation based on images from the Cassini orbiter. These authors inferred the amplitude of rotational librational motion in longitude at the orbital period by fitting a shape model to the recent Cassini ISS images. By a quasiperiodic approximation of the orbital motion, we describe how the orbital swap impacts the rotation of the satellites. To that purpose, we have developed a formalism based on quasi-periodic series with long and short-period librations. In this framework, the amplitude of the libration at the orbital period is found proportional to a term accounting for the orbital swap. We checked the analytical quasi-periodic development by performing a numerical simulation and find both results in good agreement. To complete this study, the results regarding the short-period librations are studied with the help of an adiabatic-like approach.

Introduction

The orbital motion of Janus and Epimetheus presents a peculiar horseshoe-shaped orbit resulting from a 1:1 orbital resonance (e.g. [START_REF] Dermott | The dynamics of tadpole and horseshoe orbits II. The coorbital satellites of saturn[END_REF][START_REF] Yoder | Theory of motion of Saturn's coorbiting satellites[END_REF][START_REF] Murray | Solar System Dynamics[END_REF]Jacobson et al. 2008 and references therein). Every four years the two satellites swap their orbits by a few tens of kilometers as a result of their close encounter. As the mass of Janus is 3.6 times greater than the mass of Epimetheus, the dynamical motion of the latter is more sensitive to the swap than the dynamical motion of Janus.

The rotational motion of the satellites depends mainly on the gravitational torque of Saturn acting on the dynamical figure of each moon. The expression of the gravitational torque is:

T = 3GM S r 3 u × [I] u (1)
with G the gravitational constant, M S the mass of Saturn, [I] the inertia tensor of the moon, r the distance between Saturn and the moon, and u the unit vector toward Saturn in the moon's reference frame. The gravitational torque T depends on the relative Saturn-moon distance, hence the swap also yields his signature on the rotational motion of the satellites. First estimates of the rotational motion of the two coorbital satellites Janus and Epimetheus have been obtained by [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF]. From images provided by the Cassini orbiter, they fitted a numerical shape model of the moons, which included the amplitude of the libration in longitude. The libration in longitude corresponds to the oscillation of the body along its equatorial plane. [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF] obtained an amplitude of 5.9 • ± 1.2 • for Epimetheus. For Janus, the uncertainty on the fit of the libration determination is too large to yield an accurate librational amplitude. However, based on their shape model, Tiscareno et al. suggested a value of 0.33 • ± 0.06 • for the amplitude of the libration in longitude. In addition, they identified an unexplained constant phase of 5.3 • ± 1 • for Janus, whereas for Epimetheus such offset is in the error bar. A recent numerical study by [START_REF] Noyelles | Theory of the rotation of Janus and Epimetheus[END_REF] explored the three-dimensional rotational motion of these satellites based on the numerical shape deduced by [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF]. Noyelles' study suggests a strong influence of the swap on the rotational motion of Janus and Epimetheus, which we propose to explore in the present study.

The orbital motion of Janus and Epimetheus appears to be very regular, at least over a timescale of several thousands of years. Thus, the trajectories of these satellites can be considered as quasi-periodic. Schematically these trajectories evolve on three different timescales. The shortest corresponds to the mean motion of the satellites with a period of about 0.7 days. The second component has a period of 8 years and is associated with the close encounters of the satellites. The long-period component is the secular variations of the satellites eccentricities and inclinations over periods of a few thousands years. Since the rotation of Janus and Epimetheus is synchronous with their orbital motion, it reflects these different timescales.

One of the goals of this paper is to understand the influence of the 8-year horseshoe motion on the rotational librations of Janus and Epimetheus. To this purpose, we develop an analytical solution of the rotation of these bodies that can simulate the main features of their spin. Then, in Section 3, we model the orbits of the two satellites through quasi-periodic expansions. The fourth section is dedicated to the description of the librational motion of each satellite. We develop three approaches to describe the librational motion in details: (1) a quasiperiodic development that highlight the fundamental frequencies involved in this problem, (2) an adiabatic invariant approach that focuses on the short-period librations, and (3) a numerical approach to reach high accuracy. Then we discuss the offset of Janus's orientation and investigate the effect of high spherical harmonics (order 3) and tidal coupling. We also discuss the influence of triaxiality on libration amplitudes.

Physical librations for Keplerian orbit

First, let us recall the librational response of a satellite in synchronous spin-orbit resonance with a fixed Keplerian orbit (constant semi-major axis a and eccentricity e). The position of the moon is determined through its relative distance with Saturn r, and its orbital longitude is defined by the draconic true longitude (angle between the body and its line of nodes) denoted by v. The mean anomaly will be denoted by ℓ, while the true anomaly is denoted by f . We neglect the effect of the obliquity, which is small [START_REF] Noyelles | Theory of the rotation of Janus and Epimetheus[END_REF], so that the orientation of the body is specified by the angle θ defined with respect to the line of nodes.

The dynamical equation governing the rotation of the moon is the angular momentum balance with the gravitational torque exerted by Saturn. If we introduce the physical libration γ by the relation θ = ℓ + ω + γ, where the angle γ represents the oscillations around the uniform synchronized rotation and ω the fixed argument of the pericenter, the angular momentum equation reads:

γ + σ 2 2 a r 3 sin 2(γ + ℓ -f ) = 0, (2)
The frequency σ is the frequency of the free libration, also called frequency of the proper libation. It is equal to σ = n 3(B -A)/C where A < B < C are the normalized moments of inertia of the satellite and n its mean motion.

For a small eccentricity e, the difference between the true and the mean anomalies is approximated at firstorder in e by f -ℓ = 2e sin ℓ. In addition, we approximate a = r, and, for γ small, the linearized equation of Eq. ( 2)

is γ + σ 2 γ = 2eσ 2 sin ℓ (3)
The librational solution is then simply

γ = A γ sin (σt + φ γ ) + 2eσ 2 σ 2 -n 2 sin ℓ (4)
where A γ , φ γ depend on the initial conditions and the right-hand side term is the forced libration. The forced libration oscillates at the mean motion frequency, and its amplitude is proportional to the ratio of the eccentricity to the difference between the square of the libration frequency σ and the square of the forced frequency n. Therefore, the amplitude of the forced libration depends on both the magnitude of the torque and the proximity of the free libration frequency to the orbital frequency.

In the following sections, we describe the orbits and investigate the impact of the horseshoe-shaped orbit on the physical librations.

3 Orbital description of Janus and Epimetheus

Osculating elliptical elements and fundamental frequencies

The co-orbital satellites Janus and Epimetheus are famously known to exchange their orbits every four years. This swap takes a short time-span, which does not exceed six months. In order to model the peculiar orbital motion of these satellites, we numerically integrate the three-body problem composed of Saturn, Janus, and Epimetheus, including the oblateness J 2 of Saturn. By using the frequency analysis developed by [START_REF] Laskar | Secular evolution of the solar system over 10 million years[END_REF][START_REF] Laskar | Frequency map analysis and quasiperiodic decomposition[END_REF] for the purpose of Celestial Mechanics studies, we express the numerical solution as a quasi-periodic function expanded in Fourier series, where each frequency is a linear combination (with integer coefficients) of six fundamental frequencies (proper frequencies) denoted by (n, ν, g J , g E , s J , s E ) (in the following, we use the subscript J for Janus and E for Epimetheus). The first of these frequencies, n, called proper mean motion, is associated to the mean orbital motion common to the two co-orbital satellites, and it is constant along satellite orbits. In the same way, the five other proper frequencies are also constant and might be considered as integrals of movement [START_REF] Laskar | Frequency map analysis and quasiperiodic decomposition[END_REF]. The second fundamental frequency ν corresponds to the libration frequency along the horseshoe orbits, while the four last ones are associated to the motion of the pericenters (g E , g J ) and of the ascending nodes (s E , s J ).

If the gravitational interactions between Janus and Epimetheus were negligible, then the precessions of the pericenters and of the nodes, which are only due to J 2 , would satisfy the relations g J = g E and s J = s E , and the eccentricities and inclinations of the two satellites would be constant (excepted for small short-period variations). However, due to the satellites mutual interactions, the proper frequencies of precessions of the orbits are slightly different, and both eccentricities and inclinations undergo large long-term variations whose periods are about 2760 and 7690 years, respectively. The long-term variations of the eccentricities will be neglected in the present work. Their impact on the rotation of the co-orbital satellites is mentioned in section 4.2.

The initial conditions of the numerical integrations come from the ephemeris Horizons [START_REF] Giorgini | JPL's On-Line Solar System Data Service[END_REF]. The model used in the ephemeris Horizons takes into account the gravitational interactions with the Sun and the other main satellites of the Saturnian system [START_REF] Jacobson | Revised Orbits of Saturn's Small Inner Satellites[END_REF]. Thus, in order to be consistent with that model, we fit the initial conditions of our integration in such a way that the fundamental frequencies n and ν are the same in both cases. The entire set of fundamental frequencies, as well as their associated periods, are displayed in Table 1. The frequencies (g J , g E ) and (s J , s E ) are very closed indicating that the mutual interactions are small.

Due to the high orbital precession rate generated by Saturn's oblateness, we define the elliptic elements by introducing Saturn's J 2 in the third Kepler's law, which yields:

n 2 a 3 = µ, with µ = G(M S + m J + m E ) 1 + 3 2 J 2 R S ā 2 (5)
where M S , R S are the mass and equatorial radius of Saturn, m J , m E the masses of Janus and Epimetheus, and ā the barycenter of semi-major axis Janus and Epimetheus (see Formula ( 7)). The elliptic elements (a, e, I, λ, ̟, Ω)1 are the elements of the unique ellipse tangent to the planetocentric velocity at the planetocentric location of the satellite, assuming that the third Kepler law is given by equation ( 5). As shown by Figure 1, this definition of elliptic elements removes the main orbital oscillations from the elliptical elements [START_REF] Greenberg | Apsidal precession of orbits about an oblate planet[END_REF]. In addition, it is well-known that changing the value of µ shifts the mean value of the semi-major axis by a quantity of the order of J 2 (R S /ā) 2 , which leads to the translation of about 600 km clearly visible in Fig. 1.

Analytical expression of the elliptic elements' variations

In this section, we detail the quasi-periodic expansion of the elliptic elements for both satellites that will be useful for the rotation study. According to classical theories (i.e [START_REF] Dermott | The dynamics of tadpole and horseshoe orbits II. The coorbital satellites of saturn[END_REF][START_REF] Yoder | Theory of motion of Saturn's coorbiting satellites[END_REF][START_REF] Namouni | Secular Interactions of Coorbiting Objects[END_REF], for more recent developments), the variations of the mean longitudes and semi-major axes of the co-orbital satellites are accurately approximated by the expressions:

λ J ≈ λ 0 + nt + ζ E λ r mod(2π) λ E ≈ λ 0 + nt -ζ J λ r mod(2π) (6) a J ≈ ā + ζ E a r a E ≈ ā -ζ J a r (7)
where

ζ J = m J /(m J + m E ) and ζ E = 1 -ζ J .
The variables a r and λ r represent the relative semi-major axis and mean longitudes of the satellites, that is: a r = a J -a E Figure 1: Variation of the semi-major axis of Janus and Epimetheus: the upper plot shows the temporal variation of the semi-major axes of the two satellites (solid line for Janus and dotted line for Epimetheus), when J 2 is included in the third Kepler's law. The semi-major axes oscillate around the mean value ā = 151436.9 km. When J 2 is not considered (bottom plot), large short-period oscillations are superimposed to the above-mentioned 8-year signal. The signals have a period of 2π/ν ≈ 8 years. In the latter case, the mean value of a is equal to 152024.4 km, which is about 600 km greater than when J 2 is taken into account.

and λ r = λ J -λ E . These relations clearly reflect the symmetries between the orbits of Janus and Epimetheus. Formula (7) implies that the barycenter of the semi-major axes, namely

ζ J a J + ζ E a E ,
√ a J + ζ E √ a E
becomes an integral of the motion. The relative semi-major axis ar remaining always very small with respect to a J and a E , the stated property holds.

time. In other words,

d dt (ζ J λ J + ζ E λ E ) ≈ n.
It turns out that the main variations of a and λ are given by the relative motion in coordinates (a r , λ r ).

Neglecting the terms of powers greater than two in eccentricities (which are very small for these satellites), the relative motion of the satellites satisfies the differential system3 :

   ȧr = 2εnā 1 -(2 -2 cos λ r ) -3/2 sin λ r λr = - 3n 2 a r ā , with ε = m J + m E M S + m J + m E (8)
The solutions of these equations can be expanded in a Fourier series as:

a r (t) = ā   p≥1 α (r) p cos pνt + ϕ (r) p   with α (r) 2p = 0 (9) λ r (t) = π + p≥1 β (r) p sin pνt + ϕ (r) p with β (r) 2p = 0 ( 10
)
where ν is the frequency of the relative motion, that is the second fundamental frequency defined in Section 3.1. It corresponds to a period of eight years for Janus and Epimetheus. The vanishing of the even coefficients in the series ( 9) and ( 10) is due to the symmetries of the system (8). Indeed, the invariance of the differential system by the transformation z -→ -z where z = (a r , λ r ), leads to the relations z(t + T /2) = -z(t), where T = 2π/ν is the period of the solution. Applying this relation to the Fourier expansion of z we obtain α (r) 2p = β (r) 2p = 0. While the differential system ( 8) is integrable, it is hard to get an analytical expansion of its solutions in the form of a Fourier series. Consequently, the coefficients α p have been inferred from the solutions of the numerical integration. A truncated expression of this expansion is given in Table 2. The comparison between the third and fourth columns, which display the coefficients α (r) p and β (r) p , respectively, emphasizes the different decreasing speed of theses sequences. While the coefficients of the relative semi-majors axis seem to decrease slowly, the sequence β (r) p converges more rapidly. Indeed, Equation ( 8) imposes the coefficients

β (r) p to be proportional to p -1 α (r) p .
In order to illustrate the convergence of the series ( 9) and ( 10) towards the solution of the equation of relative motion (8), we consider different approximations of these series for which the N first terms are summed. To the relative orbit of the two co-orbital satellites in the plan (a r , λ r ) visible in Fig. 2 (bold curve), is superimposed the approximated orbit (dashed curves) obtained by varying the integer N . The ellipse obtained for N = 1 provides a very crude approximation of the relative orbit, while the approximation generated at N = 30 starts matching accurately the numerical solution. In addition to the central symmetry z -→ -z mentioned above, Fig. 2 emphasizes a second symmetry with respect to the axis of coordinates. These symmetries impose relationships between the phases ϕ p , whose description is beyond the scope of this paper. The solid red curve plots the orbit deduced from the numerical simulations, while the dotted curves stands for several approximation given by formulas ( 9) and ( 10). See text for more details.

Finally, from the periodic representations ( 9) and (10) of the quantity z(t) = (a r (t), λ r (t)) and according to for-mulas ( 6) and ( 7), the trigonometric approximations of the semi-major axis and mean longitude of the satellites read:

a (x) (t) = ā   1 + 1≤p≤N α (x) p cos pνt + ϕ (x) p   = ā 1 + A (x) N (11) λ (x) (t) = λ 0 + nt + 1≤p≤N β (x) p sin pνt + ϕ (x) p = λ 0 + nt + B (x) N (12) α (x) p = (1 -ζ x )α (r) p , β (x) p = (1 -ζ x )β (r) p , ϕ (J) p = ϕ (r) p , ϕ (E) p = ϕ (r) p + π (13)
where the index x replaces J and E, whether we consider Janus or Epimetheus. The numerical values of the coefficients α

(x) p and β (x) p
are reported in the fifth and sixth columns of Table 2 for Janus and in the seventh and eighth columns for Epimetheus. In addition, the ratios of each α 2 are close to 3.6, in agreement with the formulae (13). The short-period oscillations do not appear in that table because they are negligible in comparison to the other parameters.

(E) p /α (J) p and β (E) p /β (J) p listed in Table
In the following section, we use this representation of the elliptical elements of the satellites, especially for the mean longitudes, to develop an elementary perturbation theory describing the rotation of Janus and Epimetheus. We consider in the next section that the eccentricities and inclinations are constant, as underlined in section 3.1 and discussed in Section 4.2. In the same way, the precession of the pericenters and nodes will be approximated by assuming uniform motion defined as:

e(t) = ē, ̟(t) = gt + ̟ 0 (14) I(t) = Ī, Ω(t) = st + Ω 0 (15)
4 Physical librations for Janus and Epimetheus

Perturbative analysis

Dynamical equations

The equation governing the physical libration is inferred from the angular momentum balance equation projected onto the equatorial plane of the body θ + 3 2 Table 2: Periodic approximation of the semi-major axes of Janus and Epimetheus. The relative quantities ar and λr according to Equations ( 9) and ( 10) are represented in the third and fourth columns. The fifth and sixth columns contain the coefficients of a (J ) and λ (J ) appearing in Equations ( 11) and ( 12), while the seventh and eighth columns correspond to the same quantities for Epimetheus. The frequencies listed in the second column are equal to pν where p is in the first column and ν is given in Table 1. The phases in the ninth column apply directly to Equation ( 9) and ( 10), and can be adapted to Equation ( 11) and ( 12) by means of Equation ( 13). The origin of time in formulas ( 9) is equal to 1949-Dec-28 00:00:00.0000 (JD 2433278.5) from the ephemeris Horizons.

B -A C GM S r 3 sin 2(θ -v) = 0 (16) p Frequency āα (r) p β (r) p āα (J) p β (J) p āα (E) p β (E) p ϕ (r) p (rad/day) (km) (rad) (km) (rad) (km) (rad) (rad) 1 -2.14714e-
where all variables have been defined in Section 2. Since the angle θ -v remains always small, the linearization of the equation ( 16) is a valid approximation to the rotation. Therefore, in this section, we will consider the linear timedependent equation:

θ +σ 2 a r 3 (θ -v) = 0, with σ 2 = 3 GM S a 3 (B -A) C (17)
This equation is not integrable because the quantities σ, a/r and v are implicit functions of time. However, by using the expansions in quasi-periodic functions of time presented in section 3 for the orbital elements, and an elementary perturbation theory, we can obtain the approximated solution of this equation. Let us first introduce the physical libration γ. This angle is defined as the oscillation of θ around the uniform motion (n -s)t, where n -s is the main frequency of the draconic mean longitude λ-Ω = ℓ+ω. Therefore the physical libration reads: γ = θ -(n -s)t -ℓ 0 -ω 0 , the angles ℓ 0 and ω 0 being the initial values of the mean anomaly and of the argument of the pericenter of the considered satellite4 . Consequently, using the relations ( 12), ( 14), ( 15), the definition of γ and v = f + ω for the true anomaly f , the angle θ -v also reads:

θ -v = γ -v + (n -s)t + ℓ 0 + ω 0 = γ -[ℓ + ω -(n -s)t -ℓ 0 -ω 0 ] -[f -ℓ] = γ -B N -2e sin ℓ (18)
where terms of order 2 and greater in eccentricity have been neglected. At this point, it is convenient to use the function y = γ -B N . Indeed, although the amplitudes of the terms contained in B N are very large, their frequencies are small and therefore the acceleration generated by B N , which is of order ν 2 , is negligible with respect to σ 2 . Under these approximations, Equation (17) becomes:

ÿ + σ 2 a r 3 y = 2eσ 2 a r 3 sin ℓ (19)
This equation has the same form as in the Keplerian case (see Equation ( 3)), but it is not periodic anymore but quasiperiodic on the time because, according to formulas ( 12) and ( 14), ℓ reads:

ℓ = λ -̟ = ℓ 0 + ňt + B N , with ň = n -g (20)
so = e iℓ0 e iňt 1≤q≤N k∈Z

e iℓ = e iℓ0 e iňt
J k (β q )e ik(qνt+ϕq) (21)
where the J k (x) are the Bessel functions (see Appendix A). Let us recall that in the previous expression, the coefficients β q vanish when q is even. Applying usual properties of the Bessel functions that are recalled in Appendix A, we deduce that, for q odd, the ratio between the coefficients J k (β q ) computed for Epimetheus and for Janus is well approximated by:

J k (β (E) q ) J k (β (J) q ) ≈ ζ J ζ E |k| ≈ 3.6 |k| (22)
For this reason, the coefficients of the expansion (21) decrease with increasing k much more rapidly for Janus than for Epimetheus. Consequently, the number of terms necessary to approximate e iℓ to a given accuracy using a truncated expression of ( 21) is different for the two moons (see Section 4.2.2). In order to simplify the following development, we present the series containing only the first term N = 1, which corresponds to approximate the relative orbit of the two bodies by the green ellipse in figure 2. This makes a crude simplification regarding the orbital motion of the moons but the generalization to N higher than 1 is then straightforward. Then we have:

e iℓ = e iℓ0 e iňt k∈Z J k (β 1 )e ik(νt+ϕ1) (23)
It turns out that, under these approximations and assuming as in ( 14) that the eccentricity is constant and denoted ē:

e sin ℓ = ēS with

S = J 0 (β 1 ) sin (ňt + ℓ 0 )+ p≥1 J p (β 1 ) (sin((ň + pν)t + pϕ 1 + ℓ 0 )+ (-1) p sin((ň -pν)t -pϕ 1 + ℓ 0 ) (24)
Let us mention that in the expression (23) the index of the summation k belongs to Z, while in (24) the summation is restricted to positive integers, p. Using (23) and expanding (a/r) 3 at first order in eccentricity, we also get the expression of:

a r 3 = 1 + 3ēC (25)
The expression of C is the same as S, where sine functions are replaced by cosine functions. The last term that we have to expand is σ 2 = σ2 (ā/a) 3 . By (11) with N = 1, we have a 3 ≈ ā3 (1 + 3A 1 ) and consequently:

σ 2 ≈ σ2 (1 -3A 1 ) with σ2 = 3 Gm ā3 B -A C ( 26
)
By substitution of the relations ( 24), ( 25) and ( 26) in the expression (19) this equation becomes:

ÿ+ σ2 (1-3A 1 )(1+3ēC)y = 2ēσ 2 (1-3A 1 )(1+3ēC)S (27)
If we split y into a sum of terms of decreasing magnitude as y = y 0 + y 1 + • • • , we obtain the following system of equations:

ÿ0 + σ2 y 0 = 2ēσ 2 S (28) ÿ1 + σ2 y 1 = 3σ 2 (A 1 -ēC) y 0 + 6ēσ 2 (ēC -A 1 ) S (29) . . .
The solutions of these equations are trigonometric series, whose frequencies are linear combinations with integer coefficients of fundamental frequencies of the satellites (n, ν, g, s) and of the frequency of the free libration σ. The solution independent of the free frequency σ is usually called "forced solution" and it is a quasi-periodic series of frequencies (n, ν, g, s).

When weak dissipation is introduced in the system, almost trajectories converge towards a quasiperiodic attractor (see [START_REF] Celletti | Measures of basins of attraction in spin-orbit dynamics[END_REF]). This quasiperiodic trajectory is very close to the forced solution described in the conservative system. Indeed, the fundamental frequencies of these two solutions are the same, only the amplitudes and phases are slightly different (see section 5.2). Consequently, it is relevant to focus on the forced solution in the conservative problem.

Forced librations

To begin with, let us associate to a quasi-periodic function f the function f such that:

if f (t) = p f p sin(v p t + φ p ) then f (t) = p σ2 σ2 -v 2 p f p sin(v p t + φ p )
With these notations, the general solution of (28) reads:

y 0 (t) = h sin(σt + ψ) + 2ē Ŝ ( 30
)
where h and ψ are arbitrary constants. In this section, we focus on the forced solution, so we put h = 0. As a consequence, It is easy to verify in (28) that the contributions of ēy 0 C and ē2 CS, denoted respectively ē y 0 C and ē2 CS, is a second-order in eccentricity, and that A 1 y 0 and ē A 1 S are of order ēα 1 . As, for Janus and Epimetheus, the coefficient α 1 is lower than ē, then the term y 1 can be neglected. Finally, the forced libration can be approximated by the expression:

γ(t) = 1≤q≤N β q sin (qνt + ϕ q ) + 2ēσ 2 J 0 (β 1 ) σ2 -ň2 sin (ňt + ℓ 0 )+ 2ēσ 2 p≥1 J p (β 1 ) sin((ň + pν)t + pϕ 1 + ℓ 0 ) σ2 -(ň + pν) 2 - (-1) p sin((ň -pν)t -pϕ 1 + ℓ 0 ) σ2 -(ň -pν) 2   (31)
Let us recall that, in order to give a simple expression of the short-period libration, we put N = 1 in the expression (21). For this reason, only β 1 appears in the short-period part of (31). In contrast, in the long-period component of γ (first term of Formula ( 31)), N is arbitrary. The librational angle γ is split in two types of terms exhibiting different behaviors. The first type corresponds to the 2π/ν-periodic terms that depend only on the coefficients β p , i.e., on the mean longitudes of the satellites. For these long-period terms, the dynamical figure has no influence.

The second type includes terms that vary rapidly (quasiperiodic with short frequencies ň± pν) and depend on the triaxiality of the body (B -A)/C through the libration proper frequency σ.

The amplitudes of the rapidly oscillating terms depend on the magnitude of the forcing, 2ēJ p (β 1 ) and on the proximity of the forcing frequency ň ± pν with the libration proper frequency σ. In the case of Janus, the proper frequency is 4.96 rad/day, which is far from the resonance, whereas for Epimetheus the proper frequency is equal to 8.52 rad/days, and its influence on the amplitude is substantial (see [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF].

By contrast to the Keplerian case, the swap results in the amplitude of the term associated to the frequency ň = n -g to be proportional to J 0 (β 1 ). This term is of the order of 1 for Janus, but it is significant in the case of Epimetheus as close to 0.22. Therefore, for both satellites, the rotation significantly departs from the Keplerian case.

Proper libration

In this section we investigate the proper libration of the moons (also called free libration) and we especially focus on the influence of a small divisor on the solution. By dissipative effect the proper libration is expected to be small and the damping time scale is short. However excitation mechanism might exist like recent impact for example. To focus on proper libration, let us remove the external forcing by imposing ē = 0 in (28) and (29). Then, the solution (30) reads y 0 (t) = h sin(σt + ψ), where the amplitude h is small but different from zero. Consequently, by substitution of y 0 in the equation ( 29), we get:

y 1 (t) = 3σ 2 α 1 h 2 sin((ν + σ)t + ϕ 1 + ψ) σ2 -(ν + σ) 2 - sin((ν -σ)t + ϕ 1 -ψ) σ2 -(ν -σ) 2 = - 3 2 σα 1 ν h sin(νt + ϕ 1 ) cos(σt + ψ) + O( ν σ ) (32)
Contrarily to the case of the forced libration, the amplitude of the term A 1 , which is equal to α 1 (see Table 2), plays a major role here. Indeed, in the present situation, that term is multiplied by the factor σ/ν, which is about 2200 for Janus and 3700 for Epimetheus. Thus it generates a second-order solution y 1 , whose size is comparable to the solution of order one y 0 . Then, using the values of σ given in Table 3 for Janus and Epimetheus, the librational responses for the two satellites are:

y J = h (sin(σt + ψ) + 0.3 sin(νt + ϕ 1 ) cos(σt + ψ)) (33) y E = h (sin(σt + ψ) + 1.87 sin(νt + ϕ 1 ) cos(σt + ψ)) (34)
The proper librations are combinations of a sine term with a constant amplitude and a cosine term with an amplitude varying at the swap frequency. For Janus, the main term is the sine component, whereas for Epimetheus, it is the cosine component.

Numerical study and quasiperiodic representation of the libration

Quasiperiodic decomposition

We numerically integrate the orbital and rotational dynamical equations ( 16) with the triaxialities listed in Table 3 in order to determine the accuracy of analytical solution. To clearly separate in the frequency analysis the frequencies (ň ± pν), which are quite close, we have to integrate the trajectories over a long time-span of about 400 years in the future and 400 years in the past. We also integrate during 8000 years to study the very slow variations of the eccentricities of the two satellites. This point will be discussed below. We focus on the forced libration but the initial conditions of such trajectory are not known. Hence, we use an iterative method based on the frequency analysis to converge towards this trajectory by removing the free libration amplitude (see [START_REF] Couetdic | Dynamical stability analysis of the HD202206 system and constraints to the planetary orbits[END_REF] Section 4.5).

As predicted by the theoretical approach stated in section 4.1.2, the libration angle γ = θnt -ℓ 0 -ω 0 can be naturally written as the sum of two components: γ = γ l + γ s , where γ l contains the long-period harmonics and γ s the short-period ones. The frequencies of these harmonics take the form: ň ± pν. The long-period oscillations whose amplitude is about 36 • for Janus and 130 • for Epimetheus, are very well described by the first summation in the expression (31), that is:

γ l (t) = q≥1, q odd β q sin (qνt + ϕ q )
where the numerical values of the coefficients β q appear in the sixth column (for Janus) and in the eighth one (for Epimetheus) of Table 2. In other words, the Fourier expansion of γ l is similar to the one governing the long-time oscillations of the mean longitude of the corresponding satellite.

In order to estimate the accuracy of the coefficients of the short-period terms obtained analytically in Section 4.1.2, we will focus on the Fourier decomposition of function γ s deduced from the numerical simulation.

According to the equation ( 19), the amplitude of the short-period component of the rotational libration γ s is proportional to the eccentricity of the satellite. As this quantity suffers from large very long-period variation (period of about 2760 years), the amplitude of γ s is slowly time dependent. The evolution of the eccentricities are represented in Fig. 3. During the 400-year numerical integration, the variations of the eccentricities are significant (according to Fig. 3: e J increases by ∼ 20%, while e E decreases by more than 40%). Consequently, the amplitudes of the short-period librations, which are proportional to the eccentricities, are modulated with the same ratio. As a consequence, we study the quantity γs (t) = γ s (t)e(0)/e(t). The multiplicative factor e(0)/e(t), ratio of the initial eccentricity by the eccentricity at the time t, imposes to the amplitude of γs (t) to be almost constant during the integration time. During the 200 first years of the simulation, the ratio e(t)/e(0) can be accurately fitted by the quadratic polynomial P , equal to P J (t) = 1 + 1.9839 × 10 -6 t -4.1467 × 10 -12 t 2 for Janus, and P E (t) = 1 -2.3486 × 10 -6 t -1.7633 × 10 -12 t 2 for Epimetheus. The time is counted in Julian days from 1949-Dec-28 00:00:00.0000 (JD 2433278.5). These polynomials approximate the considered ratio with a relative accuracy of 0.32% for Janus and 0.45% for Epimetheus. Consequently, the value γ s (t) of the short-period component libration at an arbitrary time t inside the 200-years considered interval can easily be deduced from γs (t) using the relation γ s (t) = P (t)γ s (t).

The main terms of the Fourier decompositions of γs are displayed in Table 4 for Janus and Table 5 for Epimetheus. The solutions are given in the form:

γ(t) = p γp sin(f p t + ψ p ) with f p = j p ň + k p ν (35)
The first and second columns contain the amplitudes of γp from the numerical simulation and the analytical expression (31). The comparison between these two columns will be discussed in the next section. The third column contains the frequencies f p , while the integers j p and k p are displayed in the fourth and fifth columns. The last column presents the phases ψ p deduced from the numerical simulation.

In the case of Janus (Table 4), the term at the orbital period is dominant, with an amplitude of 0.0052 radians, that is 0.3 • . By adding the following terms the amplitude is slightly modified. We deduce that, from the polynomial interpolation of e J given above, the amplitude of γ s increases for 0.34 • in 1950 to 0.36 • in 2010, which is very close to the value of 0.33 • given in [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF].

For Epimetheus (Table 5), no dominant term appears clearly and the decrease of the coefficients in slower than for Janus. The amplitude of the short-period oscillations given by the sum of this 13 terms deceases from 8.77 • in 1950 down to 8.35 • in 2010, which is here again, comparable the 8.9 • given in [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF] (Table 7).

Let us now briefly discuss the accuracy of the quasiperiodic approximations of the numerical solutions. Due to windowing used in the frequency analysis (see [START_REF] Laskar | Frequency map analysis and quasiperiodic decomposition[END_REF]), the method does not provide a uniform approximation of studied signal. The accuracy is generally lower on the margin of the considered interval of time. Consequently, the time span of the analysis has been chosen such that the best accuracy is obtained in the interval I = [0 : 200] years. On the interval I, the accuracy is measured with the help of the relative error ǫ p defined as:

ǫ p = Max t∈I |γ s (t) -γp s (t)| Max t∈I |γ s (t)| (36)
In his expression, γp s is the quasiperiodic approximation of γs containing the p dominant terms of the decompo-sition. As mentioned above, the value of Max|γ s (t)| is about 0.34 • for Janus and 8.77 • for Epimetheus.

The number of terms given in Tables 4 and5 is such that the relative error ǫ p is better than 10% for both satellites (Janus 8.7%, Epimetheus 9.3%), which gives an absolute error of about 0.8 • for Epimetheus and less that 0.03 • for Janus (The number of selected terms is 13 for Epimetheus and 7 for Janus). The increase in the accuracy with respect to the number of terms included in the quasiperiodic approximation is presented in Fig. 4 where log 10 (ǫ p ) is plotted against p. As it was already mentioned, the convergence of the Fourier approximation is much rapid for Janus than for Epimetheus. This figure also shows that 45 terms for Janus and 57 for Epimetheus are necessary to reduce the relative error down to 1.2%.

Finally, Fig. 5 shows the instantaneous relative error (γ s (t)γp s (t)) /Max|γ s (t)| for p = 13 (red dots) and p = 57 (green dots) in the case of Epimetheus. The results are similar for Janus. For p = 13, the discrepancy is mainly due to an over estimate of the influence of the orbital swap on the rotation (one peak every four years). The addition of the harmonics of frequency ν ± pň for high values of p tends to erase this effect.

Comparison of the analytical and numerical QP representations

The analytical and numerical short-period component of the forced solution are listed in Table 4 for Janus. The accuracy of the short-period terms is a function of the considered harmonic, whose frequency is ň ± pν, where p is an arbitrary positive integer. For p = 0, the discrepancy is about 3% and increases for increasing p, e.g. 14% for p = 2 and to 85% for p = 3. This lack of accuracy can be ascribed to neglecting the terms β q for q strictly greater than one in the expansion (21), and consequently in the analytical solutions of the rotation (31). Indeed, as it can be shown by a straightforward calculation, the amplitudes of the terms that have been neglected are given by 2J 3 (β 1 )J 1 (β 3 ) ≈ 2 × 10 -4 for p = 0, J 2 (β 1 )J 1 (β 3 ) ≈ 4 × 10 -3 for p = 1 and J 1 (β 1 )J 1 (β 3 ) ≈ 2 × 10 -1 for p = 2. These numerical values, which are deduced from Table 6, are in good agreement with the level of accuracy mentioned above.

The librational behavior of Epimetheus is reported in Table 5. Contrarily to Janus, the amplitudes of the short-period terms at ň ± ν and ň ± 2ν are greater than the terms at ň. This is due to the fact that J 1 (β 1 ) > J 2 (β 1 ) > J 0 (β 1 ) for Epimetheus, while the relation J 0 (β 1 ) > J 1 (β 1 ) > J 2 (β 1 ) holds for Janus (see Tables 6 and7). From comparing the first two columns of Table 5 we find that the accuracy on the short-period terms ob- tained for Epimetheus is worse than in the case of Janus. It reaches about 9% for p = 1 and increases to 28% for p = 4 and even more than for 80% for p ≥ 5. As for Janus, the accuracy of the analytical solution would be strongly increased if the terms related to β 3 and possibly to β 5 were taken into account.

Adiabatic-like solution of the equation of the libration in longitude

We have seen in section 4.1.1, neglecting the quadratic terms in eccentricity, that the short-period component of γ can be approximated by Eq. ( 28), which is the equation ÿ + σ2 y = 2ēσ 2 sin ℓ (37)

when N , in the series S, tends towards the infinity. To solve (28) we have expanded y in Fourier series, here we use the fact that two time scales coexist in the problem in order to solve the equation ( 37). The frequency ν being very small with respect to ň, it is possible to expend the solution in Taylor series of the small parameter ε = ν/ň. First, let us remember that, according to (20), the mean anomaly ℓ can be written: ℓ = ňt + B(νt) were B contains the long-time variations associated to the orbital swap 5 .

If we now introduce the three angles (θ 1 , θ 2 , θ 3 ), defined by θ 1 = ňt, θ 2 = νt and θ 3 = σt, we can see any solutions of (37), as a quasiperiodic function of the three previously defined angle. More precisely, this equation being linear, its solution reads:

y(θ 1 , θ 2 , θ 3 ) = ỹ(θ 1 , θ 2 ) + ŷ(θ 3 ) (38)
ŷ being periodic and ỹ quasiperiodic. The forced solution is consequently ỹ. Therefore, in order to obtain this particular solution, it is enough to restrict our analysis to the set of the quasiperiodic functions of the two variables (θ 1 , θ 2 ). Consequently, if y is an element of this set, its time derivative reads:

dy dt = ň ∂y ∂θ 1 + ν ∂y ∂θ 2 (39)
It turns out that the forced solution of the equation ( 37) is also solution of the linear partial differential equation:

ň2 ∂ 2 y ∂θ 2 1 + 2ňν ∂ 2 y ∂θ 1 ∂θ 2 +ν 2 ∂ 2 y ∂θ 2 2 + σ2 y = 2σ 2 ē sin(θ 1 + B(θ 2 )) (40)
If we now expand the solution y in Taylor series with respect to ε, that is: y = y 0 + εy 1 + • • • + ε p y p + • • • , the functions y p satisfy the sequence of differential equations:

y 0 : ň2 ∂ 2 y 0 ∂θ 2 1 + σ2 y 0 = 2σ 2 ē sin(θ 1 + B(θ 2 )) (41) y 1 : ň2 ∂ 2 y 1 ∂θ 2 1 + σ2 y 1 = -2ň 2 ∂ 2 y 0 ∂θ 1 ∂θ 2 (42) y p : ň2 ∂ 2 y p ∂θ 2 1 + σ2 y p = -2ň 2 ∂ 2 y p-1 ∂θ 1 ∂θ 2 (43) -ň2 ∂ 2 y p-2 ∂θ 2 2 , ∀p ≥ 2
Solving iteratively these equations, it is easy to show 5 The subscript N has been removed from B to indicate that the summation can be infinite that the forced solution of the libation equation reads:

y(t) = 2ēσ 2 σ2 -ň2 (1 + εU (νt, ε)) sin(ňt + B(νt)) ε 2 V (νt, ε) cos (ňt + B(νt)) (44)
where U and V are periodic functions in νt that can be expanded in power series of ε as:

U = p≥0 ε p U p+1 , and V = p≥0 ε p V p+2 (45)
The coefficients U p and V p , which depend on B and on it first p derivatives, are deduced form ( 41) to ( 43) by induction. The approximation of zero order (when ν is neglected with respect to ň):

y 0 (t) = 2ēσ 2 σ2 -ň2 sin (ňt + B(νt)) = 2ēσ 2 σ2 -ň2 sin ℓ (46)
is the solution that we intuitively get freezing the longterm temporal variations. This solution is close to the one used in [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF] for the libration, that is: A sin(n(t) t + φ), where both amplitude A and phase φ have been fitted to the observations. In this expression, the instantaneous mean motion n(t) reflects, like B(νt) in formula ( 46), the long-term variations of the "instantaneous" mean motion due to the orbital swap.

Even if the solution y 0 given by ( 46) resembles the libration in the Keplerian case (section 2), their behaviors are very different. Indeed, in the Keplerian case the mean anomaly ℓ increases linearly with the time, while in the more realistic case that we are considering, a large 8-years periodic motion in superimposes to this linear evolution (about 36 • for Janus and more than 130 • for Epimetheus). Forget this long-term variation could lead to an error on the amplitude of the libration reaching at worst 50% for Janus and 100% for Epimetheus.

In order to evaluate the accuracy of the "frozen" solution y 0 given by ( 46), we have to estimate, at least, the size of the term U 1 which appears in the expression of y 1 . Substituting y 0 and solving the differential equation ( 42), we have:

y(t) = 2ēσ 2 σ2 -ň2 1 + 2ε ň2 σ2 -ň2 B ′ (νt) × sin (ňt + B(νt)) + O(ε 2 ) (47)
where B ′ is the first derivative of B. Deducing from numerical simulations that the upper bounds of |B ′ | are respectively 1.5 for Epimetheus and 0.45 for Janus, it turns out that the addition of y 1 in the forced solution, modified the amplitudes of the "frozen" solutions of about 0.64% for Epimetheus and 0.03% for Janus, which is far lower than the accuracy of the observations.

Discussion

5.1 Higher harmonics [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF] found in the shape fitting residuals an unexplained offset in the direction of the longest figure axis of the moons, 5.2 • in the case of Janus. For Epimetheus, the detected offset is within the error bars. These authors suggested that the departure observed for Janus is due to large lateral density anomalies. We investigate such a hypothesis by assuming that the satellites shapes depart from triaxial, hydrostatic shapes, due to mass anomalies expressed at the third degree of spherical harmonics (as for the Moon e.g. [START_REF] Eckhardt | Theory of the libration of the moon[END_REF]. In this case, in addition to periodic terms and small shifts in the proper frequency, a constant term appears in the dynamical Equation leading to a constant offset of the form

γ hs3 = ζ σ2 (-15S 33 + 0.5S 31 ) (48)
where ζ = 3n 2 C R a = 0.308 rad/days 2 for Janus and 0.198 rad/days 2 for Epimetheus, where S 33 and S 31 are spherical harmonics of order 3 and degree 3 and 1, respectively. In the expression (48), we have assumed that the spherical harmonics C 13 and C 33 are negligible with respect to C 22 . Thus, the offset can be expressed as

γ J hs3 = 0.72 • (-15S 33 + 0.5S 31 ) (49)
for Janus and

γ E hs3 = 0.16 • (-15S 33 + 0.5S 31 ) (50)
for Epimetheus. As a consequence, an offset of 5.2 • requires that the combination (-15S 33 +0.5S 31 ) is about 7 for Janus, which seems very large even for rubble-piles. In this case, the shape would be far from an ellipsoid. Nevertheless, the contribution of density anomalies cannot be completely ruled out, and more complex shape models remain to be developed in order to better assess the influence of nonhydrostatic anomalies on rotation.

Tidal dissipation

A second possible origin of offset determined by [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF] might be related to the tidal torque. Indeed, Saturn raises a tidal bulge on each moon that is shifted from the planet-satellite direction due to the inelastic response of the moon material. Saturn exerts a gravitational torque on this tidal bulge and the body responds by displacing its permanent bulge so that it cancels the average saturnian torque acting on the tidal bulge. Such a displacement has been measured for the Moon [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]) and estimated for Enceladus [START_REF] Rambaux | The librational response of Enceladus[END_REF]. First, we analytically introduced the main tidal deformation as a supplementary term in equation ( 37), describing the adiabatic behavior of the moons, of the form

T = -k 2 R 5 3GM 2 S a 6 (U 11 U * 12 -U 12 U * 11 ), ( 51
)
where R is the satellite's radius, k 2 the tidal Love number, and U ij = a r 3 u i u j . The u i are the direction cosines between Saturn and the satellite. The symbol star means the position of Saturn at the constant time delay δt resulting from the dissipation. By taking advantage of the spin-orbit synchronous resonance the quantities u 1 and u * 1 are of the order of the unity and u 2 and u * 2 are small. Then, we develop u * 2 in Taylor series with respect to δt. Finally, the deformation included in the adiabatic equation Eq. ( 41) reads ÿ + 2λ ẏ + σ2 y = 2ēσ 2 sin ℓ + 4ēňλ cos ℓ (52)

with the dissipative rate

2λ = 3k 2 R 3 C ň4 Gm δt ( 53
)
where m is the mass of the satellite. We assume that the delay is constant and equal to: δt = ň-1 Q -1 . Following the perturbative analysis performed in section 4.1.1, we solve this equation by written the particular solution in the form y = y s sin ℓ + y c cos ℓ where y s is the amplitude of the in-phase term and y c is the amplitude of the out-of-phase term raises by dissipation. The last term induces a small displacement each time that ℓ = 0 mod(2π), of

y c = δ 1 k2 Q 1 + δ 2 k2 Q 2 (54) with δ 1 = 6ē C ň2 R 3 Gm ň2 σ2 -ň2 2 (55) and δ 2 = 9 C 2 ň2 R 3 Gm 2 ň2 σ2 -ň2 2 (56)
The coefficients δ 1 and δ 2 are equal to -0.011 and 0.3205 radians for Janus and -0.603 and 11.863 radians for Epimetheus. y c depends on the rheology of the satellite through k 2 /Q. This ratio depends on the internal structure and therefore on the origin of the bodies. [START_REF] Charnoz | The recent formation of Saturn's moonlets from viscous spreading of the main rings[END_REF] have suggested that many small satellites of Saturn, and especially Janus and Epimetheus, come from the accretion of ring material in the form of lumps that separate from the rings. We expect this accretion scenario to yield homogeneous bodies whose composition is mostly water ice. For the two small moons the Love number k 2 may be computed from the following relationship (e.g., McDonald 1964):

k 2 = 3 2 1 1 + 19µr 2ρgR , ( 57
)
where ρ is the mean density for the satellite, R its mean radius, g the average surface gravity, and µ r the effective shear modulus characteristic of the rubble material.

We infer from [START_REF] Goldreich | Tidal Evolution of Rubble Piles[END_REF] that the effective shear modulus of Janus is about 0.3 GPa and that of Epimetheus is 0.2 GPa. Such a low effective modulus implies an increased value of k 2 with respect to that expected for a monolith, of the order of 4.8 × 10 -4 and 3.3 × 10 -4 for Janus and Epimetheus, respectively. With regard to the dissipation factor, its value for aggregates is poorly constrained and depends on the dissipative mechanism acting inside the satellites. Then a dissipation factor of the order of 10 to 10 3 is a possible range of values for Janus and Epimetheus. The smallest dissipation factor (Q = 10) leading to the largest displacements of the axis of figure of both satellites and implies an y c about 3.02 × 10 -5 degrees for Janus and 1.14 × 10 -3 degrees for Epimetheus, which is too small to explain the offset observed by [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF].

Influence of the triaxiality on the librational amplitudes

The short-period librations are of geophysical interest because their amplitudes depend on triaxiality, as shown in Eq. ( 31). We range the triaxiality inside the error bars provided by [START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF] for Janus 0.100±0.012 and for Epimetheus [0.269 : 0.315]. For Janus, the libration amplitude depends linearly on triaxiality, while for Epimetheus this dependence is hyperbolic because the proper period and the orbital period are close to each other (see Fig. 6).

The model suggested by [START_REF] Charnoz | Physical collisions of moonlets and clumps with the Saturn's F-ring core[END_REF] leads to homogeneous bodies. Alternatively, [START_REF] Porco | Saturns Small Inner Satellites: Clues to Their Origins[END_REF] suggested that both satellites could have accreted rings particles around a core of satellite material with a lower porosity. This idea is supported by the fact that accretion models can account for the very oblate shapes of the satellites. That model implies a contrast in density between the core and a very porous outer layer. Porosity in that layer could be as large as 60 or 70%, as has been suggested for comets. There is little constraint on the thickness of that layer. Assuming Janus core is made up of solid, pure water ice, it would have a mean radius of about 70 km.

Assuming an end-member model with a solid core of water ice and a 60% porous outer layer, decreases the mean moment of inertia by about 15% with respect to the value for a homogeneous body. Unfortunately, the size of the error bars on the triaxiality from Tiscareno et al. ( 2009) is of the order of 15%, which prevents further investigation of a possible relationship between libration amplitude and internal structure.

Conclusion

In this paper we have investigated the librational motion of the co-orbital satellites Janus and Epimetheus by using three methods: (1) a perturbative technique based on quasi-periodic expansions, (2) an adiabatic invariant approach by expanding in power series of the small parameter ν/n, and (3) a numerical integration. With the perturbative technique, we have detailed the librational behavior. For both satellites the solutions are composed of long-period librations linked to the orbital swap and short-period librations related to the orbital period. We found that the amplitudes of the short-period librations depend on the magnitude of the forcing and the proximity to the resonance, as for the librations analyzed in a Keplerian framework, but also on Bessel functions of the amplitude of orbital libration of the moons mean longitudes along their horseshoe orbit. These amplitudes bear the signature of the mass distributions in the satellites and are crucial to investigate the internal structure signature of these objects. On the other hand, the amplitudes related to the long-period librations do not contain any information on the distribution of mass. The numerical integration allows us to assess the accuracy of the perturbative development. The accuracy of the analytical solution is good for Janus but poor for Epimetheus because its mass is smaller than Janus and therefore its dynamics is more perturbed. In addition, in order to obtain a compact analytical solution easy to manipulate, we have developed an adiabatic approach yielding directly the amplitude of the short periods. The analytical approaches have been developed in the most general formalism and may be applied for co-orbitals like Telesto, Calypso, Helene and Polydeuces. The adiabatic approach seems a convenient approach to fit the short-period librations to the observations.

A Expansion in Bessel functions

The Bessel functions can be defined as Fourier's coefficients of the 2π-periodic function u -→ e ix sin u where x is a real parameter, that is:

e ix sin u = +∞ k=-∞ J k (x)e iku

Figure 2 :

 2 Figure2: Relative orbit of the satellites: the X-axis represents the relative semi-major axis ar = a J -a E in km, while the Y-axis corresponds to the relative mean longitude λr = λ J -λ E in radians. The solid red curve plots the orbit deduced from the numerical simulations, while the dotted curves stands for several approximation given by formulas (9) and (10). See text for more details.

1≤q≤Ne

  iβq sin(qνt+ϕq)

Figure 3 :

 3 Figure 3: Secular variations of the eccentricities of the two satellites (solid curve for Janus, dashed curve for Epimetheus) due to their mutual gravitational interactions. The grey rectangle represents the time interval used for the numerical simulation of the rotation.

Figure 4 :

 4 Figure4: Relative error ǫp of the Fourier approximations of γ against the number of terms included in the decomposition. The " + " correspond to Janus and the " × " to Epimetheus.

  Figure 5: Relative error of the Fourier approximation of γ in the interval [0 : 200] years = [JD 2433278.5 : JD 2433278.5] for Epimetheus. The instantaneous relative error γs(t)γp s (t) /Max|γs(t)| is plotted in red for p = 13 and in green for p = 57.

Figure 6 :

 6 Figure6: Influence of the triaxiality on the amplitude of the mode n -g k (red curve) where k is equal to J for Janus (a) and E for Epimetheus (b), n -g k + ν (green curve), and n -g k -ν (blue) from Eq. (31). The black bold curves correspond to the amplitude deduced from the adiabatic invariant model (Eq. 46).

  u -ku)du

Table 1 :

 1 Fundamental frequencies characteristic of the Saturn-Janus-Epimetheus system. These frequencies are derived from a 400 years long numerical integration of the three-body problem including Saturn's oblateness.

		Freq. (rad/day) Per.
	n	9.045924658 0.69459 (days)
	ν	-0.002147139 8.01179 (yrs)
	g J	0.034948139 0.49223 (yrs)
	g E	0.034952723 0.49216 (yrs)
	s J	-0.034811959 0.49415 (yrs)
	s E	-0.034814589 0.49412 (yrs)

Table 3 :

 3 Physical properties of Janus and Epimetheus. (1)[START_REF] Tiscareno | The rotation of Janus and Epimetheus[END_REF].

		Janus	Epimetheus
	Mean Radius (1) (km)	89.5 ± 1.4	58.1 ± 1.8
	Long axis (1) (km)	101.5 ± 1.9	64.9 ± 2.0
	Intermediate axis (1) (km)	92.5 ± 1.2	57.0 ± 3.7
	Small axis (1) (km)	76.3 ± 1.2	53.1 ± 0.7
	Triaxiality (1) (B -A)/C 0.100 ± 0.012 0.296 +0.019 -0.027
	σ (rad/day)	4.95508166	8.52504483

Table 4 :

 4 Frequency analysis in the short-period librational motion of Janus. This table displays the 7 first terms of the quasiperiodic decomposition of γ(t) = e(0)/e(t)γs (t). See the text for more details.

	Amp num	Amp ana	Freq	ňJ	ν phase num
	(rad)	(rad)	(rad/days)			(rad)
	0.0052021 0.0053573	9.0109765	1	0	4.705060
	0.0014761 0.0015412	9.0131237	1	-1	3.151328
	0.0014756 0.0015433	9.0088293	1	1	3.117592
	0.0002567 0.0002159	9.0152709	1	-2	1.596407
	0.0002521 0.0002165	9.0066821	1	2	1.532845
	0.0001427 0.0000201	9.0045349	1	3	3.087729
	0.0001412 0.0000200	9.0174181	1	-3	3.183252

Table 5 :

 5 Frequency analysis in the short-period librational motion of Epimetheus. This table displays the 13 first terms of the quasiperiodic decomposition of γ(t) = e(0)/e(t)γs (t). See the text for more details.

	Amp num	Amp ana	Freq	ňE	ν phase num
	(rad)	(rad)	(rad/days)			(rad)
	0.0816497 0.0892555	9.0088248	1	1	1.356172
	0.0813811 0.0884486	9.0131192	1	-1	1.389845
	0.0618972 0.0536462	9.0152664	1	-2	2.977836
	0.0617296 0.0546295	9.0066776	1	2	2.908991
	0.0405482 0.0349633	9.0109720	1	0	6.085400
	0.0172653 0.0194522	9.0174136	1	-3	4.566668
	0.0168758 0.0199895	9.0045305	1	3	4.462410
	0.0073675 0.0052782	9.0023833	1	4	2.877236
	0.0072106 0.0050899	9.0195608	1	-4	3.006457
	0.0052737 0.0010462	9.0217080	1	-5	4.596183
	0.0051221 0.0010948	9.0002361	1	5	4.429255
	0.0024037 0.0001873	8.9980889	1	6	2.845404
	0.0023006 0.0001774	9.0238552	1	-6	3.038194

The index "J" or "E" is added to specify that the elements are related to Janus or Epimetheus if necessary.

Equivalent formulations can be found in[START_REF] Salo | The dynamics of coorbital satellite systems[END_REF] or[START_REF] Renner | Stationary Configurations for Co-orbital Satellites with Small Arbitrary Masses[END_REF]

The definition of γ is different from the one used in the Keplerian case for in the present situation the angle ℓ + ω is not proportional to the time.

Acknowledgments

The authors thank M. Tiscareno, one of the reviewers of the present paper, for his numerous comments and useful suggestions. We would also thank J. Laskar and S. Boatto for fruitful discussions. Part of this work has been conducted at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. All rights reserved. Government sponsorship acknowledged.

These functions satisfy the two following relations that we use in Section 4.1.1 : for all p ∈ N, J -p (x) = (-1) p J p (x) for all p ∈ N, J p (x) =

x p 2 p p! (1 + O(x 2 ))

In addition, we show in Tables 6 and7 the values of the main coefficients J q (β p ) that are greater than 10 -6 . These tables are useful to evaluate the accuracy of our analytical solution (see Section 4.2.2), and also to identify terms capable of increasing the accuracy of the solution.


Bruno Pinaud
email: bruno.pinaud@labri.fr

Pascale Kuntz
email: pascale.kuntz@univ-nantes.fr

GVSR: an On-Line Guide for Choosing a Graph Visualization Software

Keywords:

It is easy to find graph visualization applications for all sorts of uses. However, choosing an appropriate application may be difficult. This poster presents a website (http://gvsr.polytech.univ-nantes.fr/) built to help users to choose a program adapted to their problems. So far, this site references eighty programs and aims at helping users both in their choices and in comparing the programs. The site is also designed as a tool repository helping the community to access and compare the available tools, and benchmark new techniques and algorithms.

Introduction

The profusion of available graph visualization applications may even confuse an expert in this field. Some programs have been developed in close partnership with the scientific community (Pajek, Cytoscape), others are purely commercial, or some are general graph manipulation and visualization software (Tulip). Generally speaking, the choice of a program well-adapted to both the data and the methodology is difficult. Some books can be used as guides [START_REF]Drawing Graphs: Methods and Models[END_REF][START_REF] Mutzel | Graph Drawing Software[END_REF], and several websites present lists of programs [3,4]. However, those websites plainly list the existing software, or make them accessible through snapshots. Consequently much effort is required to compare the various programs before choosing the best one for the problem considered.

Those observations led us to develop GVSR (http://gvsr.polytech.univnantes.fr). Its added value is to offer users query about existing software based on commonly used criteria such as scalability, implementation issues or type of uses. Our objectives are to facilitate the users' choices and to compare programs with common criteria. The website also presents the programs with a uniform text-based description. This site keeps evolving and so far contains eighty various software descriptions. In addition, the site allows users to propose new programs by simply completing an enclosed form. The site is also designed as a tool repository helping the user to access and compare the available Graph Visualization tools, and benchmark new techniques and algorithms. The whole community can benefit from the ability to reproduce published results, and from comprehensive comparisons with previous work. Thus, GVSR can be seen as a contribution to improving both the accessibility and quality of graph visualization tools.

How to use the site?

The site proposes four ways to find a software : 1) a tag cloud with the software names (the more a software page is accessed, the bigger its name is displayed); 2) the "Software List" link gives access to a simple alphabetical order list; 3) the "Advanced Search" link gives access to a search engine on the software database; 4) the "Start Browsing Now" button to start navigating in a taxonomy covering all the criteria used to described software, providing a structured exploration mode of the repository.

Each software description (Fig. 1) page is made of a screenshot, general information (e.g. website, . . . ), specific information on the visualization, technical information (e.g. license(s), . . . ) and references. At the bottom of the page, one can write a comment and score the software. After validation, these information will be added on the page. 3 Future Works GVSR keeps evolving by a regular watch on the database, the addition of new programs and functionalities. We are working on an interactive visualization of the taxonomy as a graph. We also plan to directly host samples datasets to benchmark programs.

Fig. 1 .

 1 Fig. 1. Example of a software description page.