Dr Enrique García-Melendo

Esteve Duran
email: duranobs@astrogea.org

Agustín Sánchez-Lavega

Ricardo Hueso

NUMERICAL MODELS OF SATURN'S LONG-LIVED ANTICYCLONES

New measurements of the dynamical properties of the long-lived Saturn's anticyclonic vortex known as "Brown Spot" (BS), discovered during the Voyager 1 and 2 flybys in 1980-81 at latitude 43.1ºN, and model simulations using the EPIC code, have allowed us to constrain the vertical wind shear and static stability in Saturn's atmosphere (vertically from pressure levels from 10 mbar to 10 bars) at this latitude. BS dynamical parameters from Voyager images include its size as derived from cloud albedo gradient (6100 km East-West times 4300 km North-South), mean tangential velocity (45 ± 11 ms -1 at 2400 km from centre) and mean vorticity (4.0 ± 1.5 x10 -5 s -1 ), lifetime > 1 year, drift velocity (5.3 ± 0.1 ms -1 ) relative to Voyager's System III rotation rate, mean meridional atmospheric wind profile at cloud level at its latitude and interactions with nearby vortices (pair orbiting and merging). An extensive set of numerical experiments have been performed to try to reproduce this single vortex properties and its observed mergers with smaller anticyclones by varying the vertical structure of the zonal wind and adjusting the static stability of the lower stratosphere and upper troposphere. Within the context of the EPIC model atmosphere, our simulations indicate that BS's drift velocity, longevity and merging behaviour are very sensitive to these two atmospheric properties. The best results at the BS latitude occur for static stability conditions that use a Brunt-Väisäla frequency constant in the upper troposphere (from 0.5 to 10 bar) above 3.2x10 -3 s -1 and suggest that the wind speed slightly decays below the visible cloud deck from ~ 0.5 to 10 bar at a rate ∂u/∂z ~ 2-6 ms -1 per scale height. Changing the vortex latitude within the band domain introduces latitude oscillations in the vortex but not a significant meridional migration. Simulated mergers always showed orbiting

Introduction

The Voyagers 1 and 2 flybys of Saturn in 1980 and 1981 revealed the existence of a variety of "spots" with different sizes and contrasts distributed in a wide range of latitudes [START_REF] Smith | Encounter with Saturn: Voyager 1 imaging science results[END_REF][START_REF] Smith | A new look at the Saturn System: the Voyager 2 images[END_REF][START_REF] Ingersoll | Structure and dynamics of Saturn's atmosphere[END_REF]. Among them, there was a conspicuous group of oval vortices that resembled in their aspect and longevity those found in Jupiter [START_REF] Sromovsky | Voyager 2 Observations of Saturn's Northern Mid-latitude Cloud Features: Morphology, Motions and Evolution[END_REF][START_REF] Ingersoll | Structure and dynamics of Saturn's atmosphere[END_REF][START_REF] Sánchez-Lavega | New Observations and studies of Saturn's long-lived North Polar Spot[END_REF], 2000). Most of them were anticyclones, as derived from their internal rotation speed, with the largest reaching East-West dimensions between 5000 and 10000 km.

Particularly relevant were the so called "brown spots" (abbreviated in what follows as BS) located at latitude ~ 43º North within a domain of anticyclonic wind shear [START_REF] Sromovsky | Voyager 2 Observations of Saturn's Northern Mid-latitude Cloud Features: Morphology, Motions and Evolution[END_REF][START_REF] Ingersoll | Structure and dynamics of Saturn's atmosphere[END_REF][START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF]. The most conspicuous among them had about 5000 by 3000 km, with anticyclonic circumferential flow of 30 ms -1 and a mean vorticity of 5x10 -5 s -1 [START_REF] Smith | A new look at the Saturn System: the Voyager 2 images[END_REF]. [START_REF] Sromovsky | Voyager 2 Observations of Saturn's Northern Mid-latitude Cloud Features: Morphology, Motions and Evolution[END_REF] tracked these features on Voyager 1 images during ~ 20 days measuring their zonal velocity with averaged value 4.4 ± 1 ms -1 relative to System III rotation frame [START_REF] Desh | Voyager measurement of the rotation period of Saturn's magnetic field[END_REF]. A first reanalysis of the motion of this vortex from Voyager 1 and 2 was performed by [START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF] extending the tracking to 1 year (including Voyager 1 and 2 data). In this paper we extend this analysis to the highest

A C C E P T E D M

A N U S C R I P T resolution Voyager 2 images. On the other hand, the latitude band where BS formed was particularly active dynamically because of the existence of numerous other smaller spots (oval and irregular in shape) that interacted, colliding and merging after mutual orbiting [START_REF] Smith | A new look at the Saturn System: the Voyager 2 images[END_REF][START_REF] Sromovsky | Voyager 2 Observations of Saturn's Northern Mid-latitude Cloud Features: Morphology, Motions and Evolution[END_REF].

No counterpart of this kind of spots and activity was reported at this epoch in the southern hemisphere, whose visibility from Voyager spacecrafts was worse. However, this hemisphere was inspected with the Hubble Space Telescope since early 1990s at lower resolution, but still enough to detect the larger and contrasted spots. From 1994 to 2004 we recorded in the southern anticyclonic domain, between latitudes 40ºS to 44ºS and symmetrically with respect to the northern latitudes, series of dark spots that remember the "brown spots", the largest with ~ 4000 km and moving with speeds between 1 and -5 ms -1 [START_REF] Sánchez-Lavega | A strong decrease in Saturn's equatorial jet at cloud level[END_REF], 2004). These dark spots were later captured by Cassini ISS images before the spacecraft orbit insertion in July 2004, showing mutual merging behaviour when two similar spots approached [START_REF] Porco | Cassini Imaging Science: Initial Results on Saturn's Atmosphere[END_REF]. [START_REF] Vasavada | Cassini imaging of Saturn: Southern Hemisphere winds and vortices[END_REF] observed them at high-resolution using Cassini images from March-September 2004. The visual aspect and size of the two vortices present in September 2004 were similar to those of BS being probably anticyclones according to the zonal wind profile at their latitude. In addition this latitude band is prolific in forming white storms, probably convective events, observed with HST [START_REF] Sánchez-Lavega | A strong decrease in Saturn's equatorial jet at cloud level[END_REF], 2004) and Cassini ISS [START_REF] Porco | Cassini Imaging Science: Initial Results on Saturn's Atmosphere[END_REF][START_REF] Vasavada | Cassini imaging of Saturn: Southern Hemisphere winds and vortices[END_REF].

In view of the knowledge of their main properties and similarities, brown spots in the northern and southern hemispheres are at present the best candidates for the study of the dynamics of anticyclones in Saturn, and their properties can be used to constraint the dynamical structure of the upper troposphere.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

6

A variety of models have been previously used to study the longevity, stability, motions and mergers of Jovian anticyclones using a two-dimensional quasigeostrophic models and three-dimensional baroclinic multilayer models (see [START_REF] Ingersoll | Structure and dynamics of Saturn's atmosphere[END_REF][START_REF] Ingersoll | Dynamics of Jupiter's Atmosphere' in[END_REF]; and for a recent review [START_REF] Vasavada | Jovian atmospheric dynamics: an update after Galileo and Cassini[END_REF]. We apply a similar threedimensional baroclinic multilayer model to simulate a Saturnian BS. This paper deals with the report of new and improved measurements of a brown spot's main properties seen at the Voyager epoch, and used as input to perform numerical simulation of its basic characteristics and of their mutual interactions (mergers). Our aim is to use BS properties to model Saturn's anticyclones and to use BS as a probe of the dynamical state of Saturn's lower stratosphere and upper troposphere in its latitude (from 10 mbar down to about 10 bar where water clouds are theoretically predicted). We use for this purpose the EPIC code, an Atmospheric General Circulation Model well tested in previous studies of the giant planets [START_REF] Dowling | The Explicit Planetary Isentropic-Coordinate (EPIC) Atmospheric Model[END_REF]. A similar exercise was previously performed in Jupiter (Morales-Juberias et al., 2003).

Properties of the Brown Spot

Our study used archived images of Saturn obtained by the spacecraft Voyager 2 from 15 to 24 August 1981. The Voyager ISS instrument (Imaging Science System) onboard the Voyager spacecraft used wide angle (WA) and narrow angle (NA) Vidicon cameras (see [START_REF] Smith | Voyager imaging experiment[END_REF], for details on the cameras and filters). Copies of most of the images were acquired on CD-ROMs from the Planetary Data System (PDS) deep archive located at the National Space Science Data Center (NSSDC). Table 1 gives the list and dates of the images we have used.

A C C E P T E D M A N U S C R I P T [Table 1]

The image selection was based on searching for the highest-resolution images of the vortex BS in image pairs showing the same area with a minimum temporal separation ~ 10 hr. Measurements were performed on green filtered images (effective wavelength of 566 nm) because they showed the clouds with higher contrast than in other filters, at a typical spatial resolution ~ 25 -30 km/vidicon line pair. The procedure to measure positions on Voyager images is fully described in [START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF] and Legarreta and Sanchez-Lavega (2005), so we do not repeat it here. In Figure 1 we present the visual aspect of BS1 identifying the main features used as tracers of the BS flow. Wind velocity measurements were obtained by tracking the position of individual cloud features in the vortex area. Wind velocities were measured relative to Voyagers' System III radio rotation period [START_REF] Desh | Voyager measurement of the rotation period of Saturn's magnetic field[END_REF], but note that the determination of the true rotation period of Saturn is still under debate (within a range of ± 7-8 minutes) in view of Ulysses and Cassini observations (Gurnett et al., 2007 and references therein). The size and number of the cloud elements used as tracers on the vortex were scarce and only 13 features were identified moving around BS. However, many others were tracked at other longitudes at the BS latitude allowing to retrieve the zonal wind profile. The problem resides in that Saturn's anticyclones show filamentary stream-like features in their interior with few single isolated spots (as compared e. g. to Jupiter), a behaviour confirmed by Cassini images [START_REF] Vasavada | Cassini imaging of Saturn: Southern Hemisphere winds and vortices[END_REF]. The track of these features was manually performed on the side-by-side (or blinking) pair of images (Figure 1).

[Figure 1] A C C E P T E D M A N U S C R I P T
The green filter (central wavelength 566 nm) shows BS cloud morphology divided in two parts: an oval dark central core (950 km in zonal length and 700 km meridional width) surrounded by a lighter oval area (enclosed by bright filaments) that define the vortex outer edge. This are has a nearly elliptical shape with semi-major axis a = 3025 ± 100 km and semi-minor axis b = 2120 ± 100 km corresponding to a flattening ε = 1-b/a = 0.29 (or alternatively to an eccentricity e = (1-b 2 /a 2 ) 1/2 = 0.71). The vortex eccentricity changed with time (compare the two images in Figure 1) although preserving always the elliptical shape. The vortex size as seen by cloud albedo contrast is coincident within a 5% uncertainty with the motions traced by the streamlines in the fluid. In order to measure the flow around BS we first identified the vortex centre by tracing the ellipse semi-major and -minor axis. We found BS centre to be at latitude 43.1º ± 0.4º N. The second step was to measure in each image the polar position (r, θ) of the flow features relative to BS centre.

To measure the tangential velocity of each tracer, we calculate the length of the arc of the ellipse described by them during the tracking, assuming they moved following an elliptical path with the given eccentricity and semi-axes (a,b). The length of this arc is given by 2

1 ( , ) s E d θ θ θ ε θ Δ = ∫ , ( 1
)
where Δs is an incomplete elliptic integral of the second kind and

1/ 2 2 2 2 2 2 (2 ) sin cos ( , ) 1 (2 )cos r E r ε ε θ θ θ ε ε ε θ ⎡ ⎤ ⎛ ⎞ - ⎢ ⎥ = + ⎜ ⎟ - - ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ , (2)
with

A C C E P T E D M A N U S C R I P T 2 2 2 2 2 2 ( ) sin cos a b r a b θ θ θ = + . (3)
The tangential velocity is simply given by V T = Δs/Δt and the integral (1) was numerically calculated for each tie-point. This method gives essentially the same result as that determined by measuring V T from the radial and orthogonal components of the velocity (polar coordinates) of the tracers. In Table 2 we give the positions and velocity of the most reliable tracers of the flow shown in figure 1. Some other features where traced just outside the vortex, for example in the stagnant points at the vertex of the semi-major axis. On the average we find that V T = 45 ± 11 ms -1 at a distance r = 2400 km. A fit to the radial dependence of the tangential velocity (assuming no polar angle dependence) gives: 2 1021 0.88 0.00018

T V r r = - + - , (4)
for r given in km and V T in ms -1 . The maximum speed from this parabolic fit is 52 ms -1 .

[Table 2]

The mean vorticity, averaged over the area enclosed, is the circulation divided by the area enclosed (Hess, 1969)

ab L V ab l d V e T T π π ς ≈ • = ∫ r r , ( 5
)
where L e is the length (perimeter) of the ellipse which can be approached using the Ramanujan formula by This gives an averaged vorticity ζ = 4.0 ± 1.5x10 -5 s -1 , which is very close to the 3.3x10 -5 s -1 averaged ambient relative vorticity of the anticyclonic domain bounded in the North by an intense eastward jet stream at 47ºN (u ~ 150 ms -1 ) and in the South by a weak westward jet at 39ºN (u ~ -20 ms -1 ) (Figure 2) [START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF]. The averaged zonal flow velocity at the latitude of BS was <u (43.1ºN)> = 20.6 ± 16.0 ms -1 [START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF] and therefore, BS moved westward relative to the ambient flow with a relative speed u BS = -15.3 ms -1 .

[Figure 2]

The EPIC model atmosphere

By working in isentropic coordinates and assuming an adiabatic atmosphere, the EPIC ("Explicit Planetary Isentropic-Coordinate") numerical model [START_REF] Dowling | The Explicit Planetary Isentropic-Coordinate (EPIC) Atmospheric Model[END_REF]) solves the fully non-linear hydrostatic equations for the horizontal velocity components. The model computes, among other variables, the potential vorticity q, which is a conservative quantity that can be used as an atmospheric tracer in an inviscid and adiabatic flow [START_REF] Pedlosky | Geophysical fluid dynamics[END_REF].

As performed in our previous simulations in Jupiter [START_REF] García-Melendo | Jupiter's 24° N highest speed jet: Vertical structure deduced from nonlinear simulations of a large-amplitude natural disturbance[END_REF], the main free parameters of the simulated atmosphere model relate to the assumed separable meridional u(φ) and vertical u(z) wind profiles, and the vertical thermal profile. Only the meridional profile of the zonal wind at cloud top level is known in the simulated horizontal domain and is taken from Sánchez-Lavega et al.

( 2000). The vertical wind shear and thermal profile are mostly unknown quantities below the visible cloud deck and are the basic properties to be retrieved.

A C C E P T E D M

A N U S C R I P T

The thermal vertical profile

The thermal profile is based on the radio occultation experiments performed by the Voyager mission [START_REF] Lindal | The atmosphere of Saturn: an analysis of the Voyager radio occultation measurements[END_REF]. Initially we extended it down to P = 10bar, assuming that in its upper part Saturn's atmosphere is a mixture of non-condensable gases H 2 , He, and CH 4 , and the condensable species NH 3 , NH 4 SH, and H 2 O, after applying the thermo-chemical models proposed by [START_REF] Atreya | A comparison of the atmospheres of Jupiter and Saturn: deep atmospheric composition, cloud structure, vertical mixing, and origin[END_REF]. The thermal profile was corrected for a helium molar fraction of 0.13 according to the re-analysis performed by [START_REF] Conrath | Saturn Helium Abundance: A Reanalysis of Voyager Measurements[END_REF] of Voyager data and by assuming that hydrogen was in ortho-para equilibrium. This profile tends towards an adiabat near the 1 bar level, although it is slightly unstable (N 2 < 0, being N the Brunt-Väisälä frequency)

for P > 400mbar (Figure 3). Nevertheless, the presence of vigorous convective phenomena and long lived vortices at mid latitudes in both hemispheres [START_REF] Sánchez-Lavega | A strong decrease in Saturn's equatorial jet at cloud level[END_REF], 2004;[START_REF] Porco | Cassini Imaging Science: Initial Results on Saturn's Atmosphere[END_REF] suggests that statically stable as well as unstable areas must coexist in these bands of Saturn's atmosphere. We adopted the hypothesis that on the average the atmosphere in this region is statically stable since long-lived vortices form at all longitudes and that it becomes marginally unstable in short periods and small regions where convective storms form.

We adopted a simple vertical thermal structure by keeping the measured N(P)

profile in the upper part of the model atmosphere and making it constant below a fixed pressure value. To make the thermal profile continuous we adopted the approach used by LeBeau and Dowling (1998) in their simulations of Neptune's Dark Spot. According to their derivation and notation, assuming uniform stratification (N 2 = N 0 2 = constant), an ideal gas behaviour, and by making use of the hydrostatic approximation, the T(P)

profile writes for P ≥ P 0

A C C E P T E D M A N U S C R I P T 0 0 0 0 1 T P P T T T P P ) P ( T N N ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = κ κ . ( 7
)
Here κ is the ratio of the gas constant to the specific heat at constant pressure and T N ≡ κg 2 / N 0 2 R gas is fixed by the chosen value of both N 0 2 and the gas mixture, and governs the slope of T(P). The pair (P 0 ,T 0 ) is selected to assure temperature continuity between the stratospheric observed T(P) profile and the constant N(P) region bellow. Although for a given T N the pair (P 0 ,T 0 ) can be arbitrarily chosen, if we want a reasonable match between our model T(P) and the extended Saturn's profile, we may think of preserving the temperature minimum (the tropopause), and avoiding an excessive departure of our model T(P) from the extended Saturn's profile. These two conditions restrict the N 0 2 values that can be used. In our simulations we tested the two values of N 0 2 = 0.03x10 -4 s - 2 and 0.1x10 -4 s -2 . The resulting T(P) and N 2 (P) profiles are depicted in Figure 3.

[Figure 3]

The wind profile

We use the Voyager meridional wind profile at cloud top level as measured by [START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF] in our EPIC simulations. BS is located in the anticyclonic, south side, 47ºN jet stream. Figure 2 depicts the part of the jet stream used in the simulations including its uncertainty due to navigation and measurement errors.

To deal with the vertical behaviour of the wind profile, we assumed that U(φ,P) was the product of separate functions of latitude φ and pressure P. This is in

A C C E P T E D M A N U S C R I P T
agreement with observed profiles in the upper troposphere and stratospheric thermal winds [START_REF] Flasar | Temperatures, Winds, and Composition in the Saturnian System[END_REF] and it constitutes a convenient approximation to the unknown lower tropospheric winds. Therefore we use

U(φ,P) = u h (φ)u v (P) (8)
where u h (φ) is the zonal wind profile at the cloud level P 0 and u v (P) is a non dimensional function normalized at the cloud level pressure P 0 (P 0 ~ 500 mbar, u v (P 0 ) = 1), which modulates the amplitude of u h (φ). As in García-Melendo et al. ( 2005), we adopted a piecewise structure for u v (P) in two linear segments (see Figure 4). For P < P 0 , and assuming an m 1 e-folding vertical scale height, the expression for u v (P) is

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = 0 1 1 P P ln m ) P ( u v . (9)
[Figure 4]

Our models were run for two different values of m 1 . [START_REF] Ingersoll | Structure and dynamics of Saturn's atmosphere[END_REF] found an average value of m 1 = 0.2 for the wind decay above the cloud deck based on the results published by [START_REF] Conrath | Thermal structure of Saturn from Voyager infrared measurements -Implications for atmospheric dynamics[END_REF]. More recent results based on Cassini observations suggest that at midlatitudes there is almost no wind decay for Saturn's atmospheric layers above the 500 mbar level [START_REF] Flasar | Temperatures, Winds, and Composition in the Saturnian System[END_REF]. Although these Cassini results are limited to the planet's southern hemisphere, we also adopted as a second test value m 1 = 0.

A C C E P T E D M A N U S C R I P T

We have no information on the vertical wind shear below the ammonia upper cloud limit. Here is where we introduced our two-piece linear segment model to test different wind shear tendencies. For P > P 0 , we extended this profile to a level P 1 where we switched to slope m 2 , a parameter which could take positive, zero, or negative values. Below P 1 (higher pressure levels), we adopt the following expression for u v (P):

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = 1 2 1 P P ln m u ) P ( u v ( 10
)
) P / P ln( m u 0 1 1 1 1 + = (11)
Therefore, m 1 , N, P 1 and m 2 are the basic atmospheric parameters to be retrieved from our simulations, although, m 1 and N are somewhat constrained by the observations.

Simulation domain

Simulations were carried out in a channel spanning 60º in longitude by 15º in latitude with a grid of 256x64 points (in some cases 128x32, which represents the minimum number of points required to reproduce the vortex behaviour) for the horizontal domain. These values assured an equal latitudinal and longitudinal resolution of 0.24º per grid point, or equivalently 235km (respectively 0.47º or 470km). The channel covered the planetographic latitude interval from 36ºN to 51ºN. Vortex pair interactions were simulated at a grid resolution of 0.24º in the same horizontal domain. (2005). In this case, and as it can be appreciated in Figure 3, the bottom layer is an abysm with almost neutral stability.

Single vortex simulations

In the first set of simulations aimed to test the effects of the vertical wind shear and thermal profile on the vortex drift velocity and lifetime, we introduced a single vortex in the simulation domain. After specifying the vortex tangential velocity, its position, and its horizontal and vertical extension, EPIC introduces a geostrophically balanced vortex by modifying the initial (t=0) Montgomery potential M=C p T+gz, which appears as the forcing term in the isentropic-coordinate non-linear equations [START_REF] Dowling | The Explicit Planetary Isentropic-Coordinate (EPIC) Atmospheric Model[END_REF]. Afterwards, it is let to evolve freely (see [START_REF] Morales-Juberías | EPIC simulations of the merger of Jupiter's White Ovals BE and FA: Altitude dependent behavior[END_REF]. The vortex was introduced at latitude 43ºN with initial tangential winds of 50 and 100 ms -1 . Although this last value for V T is almost as double as the one measured from Voyager images, we wanted to make sure that the simulated vortex, which slowly damps out as the simulation goes on, acquired tangential velocities similar to the measured V T by the time it reached half of its lifetime span. In this way we also tested different values for the "hyperviscosity parameter" (see Table 3 andsection characteristic horizontal frequency [START_REF] Conrath | Termal structure and dynamics of the Jovian atmosphere 2. Visible cloud features[END_REF]. For the different possible values of N and the scale height H, this estimated ratio may take any value between 2.0 and 4.5 scale heights, so we made simulations using vertical extents going upward from 2 to 3.5 and downwards from 3.0 to 4.0 as a possible working scale heights above and below BS visible clouds.

Table 3 summarizes the range of values used for each parameter in the simulations.

[Table 3]

Results

The main goal of this set of experiments was to identify which set of values of the parameter space resulted in a vortex drift rate of ~ 5 ms -1 with respect to System III, and if it was also possible to maintain stable vortices with longevities longer than 1 year when the parameter combinations yielded good drift rates. For the physical model, we tested all the combinations for P 0 , P 1 , m 1 , m 2 , and N(z), in addition to other parameter combinations directly related to the numerical computations as listed in Table 3.

Drift velocities and vortex longevity

Vortex drift velocities turned out to be very sensitive to the different combinations of m 1 , m 2 , and P 1 . Our simulations suggested that drift rates could serve as a discriminator among the different model atmospheres. atmospheres, implied at most a +/-0.1 uncertainty for the m 2 parameter. Therefore, wind profile errors yielded drift velocity alterations below our numerical model resolution.

A second source of uncertainties is the effect of the grid resolution and hyperviscosity on drift rate and vortex longevity. We found that the drift velocity, for any particular model atmosphere, was almost no sensitive to these parameters. Hyperviscosity, a parameter introduced to avoid growing numerical instabilities, can damp vortices out too quickly if it has a too high value. It is also directly related to grid resolution, and it can easily be decreased with increasing grid resolution. High grid resolution not only represents with much better detail any simulated phenomenon, but also allows using lower hyperviscosity values. The price to pay is the need for enormously longer CPU time. For instance, doubling grid resolution from 128x32 to 256x64, means the need of almost an order of magnitude longer CPU time, which is absolutely prohibitive if we want to fully explore all the combinations of the model atmosphere parameters. We performed some high and low resolution simulations with the same atmospheric parameters, and look at the impact on drift velocity. If drift velocity was not affected then we would explore all the parameter space at lower resolution. That was the situation, drift velocities derived from high and low resolution simulation were virtually indistinguishable for the same model atmosphere, and it was not affected for different values of hyperviscosity either.

Finally, vortex drift velocity was also almost independent of the simulation time, but vortex longevity was affected by grid resolution. Vortices had longer lifetimes for high resolution and low hyperviscosity values. A way of measuring the vortex longevity is simply looking at the potential vorticity field to see when the vortex dissolves. Another

A C C E P T E D M A N U S C R I P T

way of measuring longevity is looking at how the vortex tangential velocity evolves with simulation time. As an example, Figure 6 shows the evolution of maximum vortex tangential velocity with simulation time for the parameter combination m 1 =0.0, m 2 =-0.2, P 1 =1000mbar, and N 2 =0.1x10 -4 s -2 for two different grid resolutions. For the low resolution case (grey line), vortex tangential velocity slows down at the beginning after vortex introduction, reaches a stable tangential velocity and afterwards, during dissipation, it starts oscillating until the vortex dissipates in the environment atmosphere. In the case of higher resolution, the tangential velocity remains stable with a very slow damping for a longer period of time. Despite the alterations of tangential velocity in low resolution simulations, the drift rates remained completely stable during the whole vortex lifetime. Low resolution simulations also provide a lower limit to the vortex longevity for each tested atmosphere. As an example, Figure 7 shows the evolution of drift velocities for all the simulated values of m 2 for the case m 1 =0.0, P 1 =1000mbar, N 2 =0.1x10 -4 s -2 and a 128x32 grid. They remained stable within 2 to 3 ms -1 until vortices dissolved, which again is well below the model parameter resolution.

[Figure 6]

[Figure 7]
Drift velocities close to the experimental value must be associated to longevities above 1 year, a lower bound imposed by observations. Figure 8 shows the measured longevities for the different combinations of the parameter space for the low resolution simulations. If we take into account that higher resolution grids result in simulated vortices with longer longevities, Figure 8 [Figure 8]

Global single vortex test results

We found model atmospheres where simulated vortices matched reasonably well the observed properties. In all cases the best results were obtained for negative values of m 2 . Although it is difficult to asses which is the "best" atmosphere model there are some trends resulting from the simulations. First, it is not possible to obtain drift velocities of 5ms -1 for m 1 = 0.2 and N 2 =0.03x10 -4 s -2 and only marginally for N 2 =0.1x10 -4 s -2 . This, of course, does not invalidate atmospheres with vertical wind shear above the visible cloud deck, but the better results for m 1 = 0 are more consistent with the Cassini results reported by [START_REF] Flasar | Temperatures, Winds, and Composition in the Saturnian System[END_REF]. Second, lower values of P 1 seem to favour vortices with a longer lifetime and good drift values. And third, vortices introduced in atmospheres with negative m 2 values show an aspect ratio between the north-south semi-axis closer to the measured values for BS, while for progressively larger positive values of m 2 , the vortices tend to be more elongated.

Latitude migration experiments

Since the meridional size of the vortex is smaller than the latitudinal width of the band comprised between the westward jet at 39ºN and the eastward jet at 47ºN, a latitude migration of BS was in principle possible, so we explore this possibility in our experiments. In the first simulation we introduced the anticyclone at φ = 45ºN and in the second at 41ºN, i.e. displaced ± 2º from its nominal latitude, using a model atmosphere

A C C E P T E D M A N U S C R I P T
with a high static stability (N 0 2 = 0.1x10 -4 s -2 ), and m 1 = 0.0, m 2 = -0.2 and P 1 = 500mbar.

Figure 9 shows that there is no evidence of a large migration for the two studied cases.

The vortex introduced at 41ºN just oscillated around this latitude. When introduced at 45ºN it migrated about one degree to the south during the initial first days, staying afterwards about 1 degree to the north of the measured BS position.

[Figure 9]

Vortex pair interactions

There is still a lack of detailed observations about the dynamical processes that occur during vortex merging in Saturn. [START_REF] Smith | A new look at the Saturn System: the Voyager 2 images[END_REF] and [START_REF] Sromovsky | Voyager 2 Observations of Saturn's Northern Mid-latitude Cloud Features: Morphology, Motions and Evolution[END_REF] showed examples of mergers between anticyclones that occurred in the latitude band of BS during the Voyager encounters. [START_REF] Porco | Cassini Imaging Science: Initial Results on Saturn's Atmosphere[END_REF] observed in Saturn's southern hemisphere using Cassini low resolution images, how two vortices spiralled toward each other during the merging process that took place in about three to four days.

In this new set of experiments we simulated the interaction between two vortices at high resolution to perform two types of tests. In the first one, two identical vortices were introduced in the domain separated by a small distance in latitude, one of them at 43ºN, the other typically one degree northwards. They were initially separated in longitude by about 40 degrees to make sure that the interaction would only occur after the short initial period of geostrophic adjustment of the EPIC model. In the second experiment, one of the vortices is introduced at 43ºN with the same parameters as in the single-vortex experiments, but the second vortex is smaller, typically with dimensions about two-thirds that of the big one. Again, they are placed with a small separation in latitude, but a large one in longitude.

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT

22

We chose the same model atmosphere used to perform the vortex migration experiments. When vortices settled after their introduction, they approached each other showing a centre-to-centre separation of 3 degrees. Figure 10 illustrates the results. The left column shows our simulations when both vortices are identical. During the interaction process, orbiting started to take place when both vortices were ~5 longitude degrees apart, they then quickly spiralled around each other in a symmetric merging to produce a single vortex in about three days, which agrees with the temporal scale observed for the southern hemisphere vortex mergers [START_REF] Porco | Cassini Imaging Science: Initial Results on Saturn's Atmosphere[END_REF]. The case of two vortices with different sizes is shown in the right column of Figure 9. In this case, the small vortex starts orbiting the big one when they are between 7 and 10 degrees apart in longitude, resulting in a mixing process that again, takes about three days. In this latter case it is interesting to note that the merging process is accompanied by the formation of potential vorticity filaments (i. e., fluid filaments) linking both spots.

These experiments strongly suggest that orbiting is a phenomenon accompanying a two vortex merger, and that the time required for the two vortices to merge as well as the orbiting phenomena depend on the flow pattern and relative speed of the two vortices.

[Figure 10]

Discussion

Our measurements clearly indicate that the BS vortex is in geostrophic balance with Rossby number Ro = V T /f L ∼ 0.05, being the planetary vorticity f = 2Ω sinϕ = 2.2x10 -4 s -1 at the BS latitude. The vortex size L is above the Rossby deformation radius L D > Nh/f ~ 2200 km when using for h a vortex vertical extent ~ 4H = 150 km and for N > 3.2x10 -3 s -1 as obtained from our simulations (see section 3). Its meridional extent nearly fills the domain between opposed jet streams (see figure 2), being constrained to

A C C E P T E D M A N U S C R I P T
the North by the "ribbon wave" [START_REF] Sromovsky | Voyager 2 Observations of Saturn's Northern Mid-latitude Cloud Features: Morphology, Motions and Evolution[END_REF]. Accordingly, the Burger number (L D /L) 2 ~ 0.08 suggests that the vortex is stable against its tendency to be broken apart vertically, i. e. the density stratification in the vertical is low enough for the flow to be dominated horizontally by the planetary rotation (see e. g. [START_REF] Ingersoll | Structure and dynamics of Saturn's atmosphere[END_REF] . This is consistent with its lifetime > 1 year that is well above its characteristic dynamic time τ dyn = L/V T and destruction time by the shear wind τ shear = 1/(∂u/∂y) ~ 1 day.

In Figure 11 we show how the zonal flow behaves in terms of a simple barotropic instability criterion at the latitude of BS. The magnitude of interest β -∂ 2 u/∂y 2 changes sign close to the BS zonal wind domain so the necessary condition (Rayleigh -Kuo criterion) for the barotropic instability to occur is met [START_REF] Holton | An Introduction to Dynamical Meteorology[END_REF].

In the more realistic situation we have simulated with EPIC, the baroclinic term should be considered to study the flow stability. However we do not consider this analysis due to the sensibility of this term to the uncertainties we have in the retrieved values for N(z) and u(z) below clouds. However we can address how large could be the contribution of the vertical wind shear to the vorticity of BS which is ~ (f/N) 2 (∂u/∂z) ~ 10 -6 -5x10 -7 s -1 (using from our simulations ∂u/∂z ~ 10 -4 s -1 and a lower limit for N > 3x10 -3 ). This is low and below the contribution of the horizontal wind shear vorticity (∂u/∂y ~ 3x10 -5 s -1 ) and that of the planetary vorticity f to the averaged BS vorticity ζ = 4x10 -5 s -1 .

Conclusions

Our suite of anticyclone simulations for the main atmospheric model parameters N 0 2 , m 1 , m 2 , and P 1 , suggests that in order to obtain stable vortices with lifetimes long enough to at least match the one-year observed survival for BS, and to get drift

A C C E P T E D M A N U S C R I P T
velocities close to 5 ms -1 , winds must slightly decay with depth down to the 10 bar level at the BS latitude. A more stable thermal profile (N 0 2 = 0.1x10 -4 s -2 below log(P) ~ -1.0) yields more stable vortices. Simulations seem also favour constant winds or nearly so (m 1 ≅ 0.0) above the visible cloud deck (P 0 ≅ 500mbar), a result which seems to be consistent with the reported Cassini measurements of the thermal winds by [START_REF] Flasar | Temperatures, Winds, and Composition in the Saturnian System[END_REF]. Finally, the value of P 1 (the level for the assumed inflexion in the vertical wind profile) is constrained by the zonal velocity results (see Figure 5), and locates it between the 500 to 1000 mbar level.

The overall picture from our model simulations at the BS latitude, and within the coarseness of the values of the parameter space, can be summarized in the following way. Saturn's upper troposphere at mid-latitudes has almost no wind shear above the 1 bar level. Below that level there is a slight wind decay with a vertical wind shear ∂u/∂z ~ 2 to 6 ms -1 per scale height in a slightly statically stable environment at the BS latitude. We would like to emphasize that our inferences on N(z) and u(z) are obtained within the context of the EPIC model and of the reference atmosphere we have used. and N 2 (see Table 3). Only decreasing winds with depth resulted in drift velocities close to the experimentally measured values (horizontal line). The dark dot in panel B marks the drift velocity for m 1 =0.0, m 2 =-0.2, P 0 =500mbar, P 1 =1000mbar, and N 2 =0.1x10 -4 s -2

for the nominal zonal wind profile. When the nominal zonal profile is substituted by its upper and lower error limits, the resulting vortex drift velocities are marked by the light dots. 3. P 0 is the cloud top level (500mbar) where u v (P 0 ) = 1.

Figure 5.

Simulated zonal vortex speeds for the different parameter combinations (P 0 = 500mbar) with a grid resolution of 0.47º. A) m 1 =0.0, N 2 =0.1x10 -4 s -2 . B) m 1 =0.0, N 2 =0.03x10 -4 s -2 .

C) m 1 =0.2, N 2 =0.1x10 -4 s -2 . D) m 1 =0.2, N 2 =0.03x10 -4 s -2 . The horizontal line marks the measured BS1's drift speed.

Figure 6.

Evolution of the maximum tangential velocity as function of the simulated time and spatial resolution for a vortex with the parameter combination m 1 =0.0, m 2 =-0.2, P 1 =1000mbar, and N 2 =0.1x10 -4 s -2 . The black and grey lines are the results for grid resolutions: 256x64 points (black line, 0.24º pix -1 ), 128x32 points (grey line 0.47º pix -1 ).

A C C E P T E D M

A N U S C R I P T Vortex drift velocity for different values of the vertical wind shear parameter m 2 in the case: m 1 =0.0, p 1 =1000mbar, and N 2 =0.1x10 -4 s -2 . This graph is representative of other combinations of m 1 , p 1 , and N 2 , and shows that the drift velocity is a stable parameter along the life of the simulated vortex.

Figure 8.

Vortex longevity for different combinations of P 1 , m 1 , m 2 , and N 0 2 (P 0 is always fixed to 500 mbar), and a grid resolution of 128x32. The legend bar in the bottom part of the figure codes the maximum lifetime of the vortex in days, black is for a longevity up to 100 days, while white represent longevities longer than 900 days. For several cases of negative m 2 , and specially for the case N 0 2 = 0.03x10 -4 s -2 , the resulting model atmosphere was unphysical, which means that isentropic surfaces crossed each other as a result of excessive vertical shear of the zonal wind to be consistent with hydrostatic and gradient-wind balance.

  The used time step length Δt was of 60 seconds for the low resolution simulations and 30 seconds in the high resolution experiments, well below the average Courant number limits of 100 and 55 seconds respectively. The vertical domain extends from P top = 10 mbar to P bottom = 10 bar and was divided in eight layers as described in García-Melendo A C C E P T E D M A N U S C R I P T et al.

A

  4.1). It wasgiven the adequate initial size parameters to mimic the actual measured dimensions during its mature stage, again taken at half lifetime. With almost no difference in the results, we introduced the vortex at two different levels P 0 = 1 bar and P 0 = 1.5 bar where the putative ammonia cloud forms. Our vortex has the form of a Gaussian perturbation with a limited vertical extent of a few vertical scale heights. It is, therefore, thin relative to its much larger horizontal dimensions. The vortex vertical extension under the cloud deck is unknown, but the ratio between the vertical and horizontal scales should be of the order of ~ f/(2N), where f is the Coriolis parameter as a

[Figure 5 ]

 5 Figure5shows that and increment, positive or negative, of 0.2 around m 2 = -0.2, introduces velocity variations of ±15ms -1 . Positive values for the vertical wind shear below P 0 seem to be completely incompatible with our observations as they turn out into too high velocities. Panel A of Figure2displays the part of Saturn's zonal wind profile used in the simulations and, as solid dots, all the simulated drift velocities for BS1 represented in Figure5.

  represents a lower bound which indicates that long-lived vortices in our model atmosphere are compatible with the adequate drift

  Figure captions

Figure 2

 2 Figure 2Panels A and B show the zonal wind profile (grey band indicates the uncertainties) as measured by[START_REF] Sánchez-Lavega | Saturn's zonal winds at cloud level[END_REF] that is used in the numerical simulations (nominal zonal wind profile). Solid dots in panel A mark the drift velocities of the simulated vortices for all the model atmosphere parameter combinations m 1 , m 2 , P 0 , P 1 ,

Figure 3 A

 3 Figure3

Figure 4

 4 Figure 4Vertical structure scheme of the zonal wind amplitude functions u v (P) used in our simulations. Two main cases are shown (grey m 1 =0.0 and solid lines m 1 =0.2), each one divided in the bottom part in the sub-cases corresponding to different values of the slope m 2 as specified in Table3. P 0 is the cloud top level (500mbar) where u v (P 0 ) = 1.

  Figure 7.

Figure 9 .

 9 Figure 9. Latitude position of two vortices, one introduced at 45ºN (upper line) and the other at 41ºN (bottom line) as a function of time for N 0 2 = 0.1x10 -4 s -2 , m 1 = 0.0, m 2 = -0.2, and P 1 = 500mbar. The grey band represents de interval of latitude positions for BS measured by Sánchez-Lavega et al. (2000).

  Figure 10.

Figure 11

 11 Figure 11 Upper panel: Vorticity of the mean zonal flow (meridional profile) as derived from the observations (Sanchez-Laverga et al. 2000). Lower panel: representation of the difference between the planetary vorticity β and ∂ 2 u/∂ 2 y as a function of latitude. The vertical grey line marks the position of BS1. The light grey ellipsoid represents the vortex BS-1 at the latitude scale.

  Figure 1

Table 1 Voyager 2 images selected to study the "Brown Spot" (BS) anticyclone

 1

	PDS Volume	Image Number	Date	Time
	4	C4381308	August 19, 1981	17:24:47
	4	C4382827	August 20, 1981	05:39:59
	4	C4383930	August 20, 1981	14:30:23
	4	C4385247	August 20, 1981	01:07:59
	4	C4386551	August 21, 1981	11:35:11
	4	C4387908	August 21, 1981	22:12:47
	5	C4383401	August 23, 1981	18:07:11
	5	C4384649	August 24, 1981	04:21:35

Table 2 Tangential velocity of the flow around the "Brown Spot" (BS) anticyclone

 2

	Feature	Position angle	Radial distance	Tangential Velocity
		<θ>	<r>	V T (ms -1 )
	A	80.2	2227	62.6
	B	47.1	2583	41.2
	C	11.3	2266	55.9
	D	330.1	2632	44
	E	311.5	2323	38
	F	185.6	2204	14.5
	G	218.9	2518	59.2
	H	323.4	2258	53.8
	I	43.9	2153	33.4
	J	0.5	2800	31.3
	K	358	2076	30.4
	L	224.4	2212	52.6
	M	345	2670	30.4

Table 3 Range of values of the parameters used in the simulations

 3

	Parameter	Range of values
	P 0 (mbar)	500
	P 1 (mbar)	500, 1000, 1500, 2500
	m 1	0.0, +0.2
	m 2	-0.4, -0.3, -0.2, -0.1, 0.0, +0.1, +0.2, +0.3
	Cup (scale heights)	2.0, 2.5, 3.5
	Cdown (scale heights)	3.0, 3.5, 4.0
	Vortex height (mbar)	1000, 1500
	Initial Vt (ms -1 )	50, 100
	Vortex size (a,b)	6.5º x 3.7º
	N 2 (z) (minimum) (10 -4 s -2 )	0.03, 0.1
	Domain	Long: (-30º a +30º), Lat (+36º a +51º)
	Number of layers	8
	x-y resolution (points)	128x32, 256x64
	ν 6 (10 27 m 6 s -1 )	0.18, 2.2, 2.8, 5.5, 8.8
	Δt (seconds)	30, 60

Acknowledgements

This work has been supported by the research projects AYA2003-03216, FEDER and Grupos UPV 15946/2004. RH acknowledges a "Ramón y Cajal" contract from MEC.

We acknowledge the National Space Data Center through the World Data Center A for Rockets and Satellites for providing the Voyager images. The Voyager Experiment Team Leader was Dr. Bradford A. Smith. This research also made use of the computing facilities at CESCA in Barcelona with the help of the Ministerio de Educación y Ciencia.


Keywords: Saturn, Satellites; Satellites, Dynamics

The aim of this work is to understand the absence of objects along the orbits of Mimas and Enceladus in contrast to their presence at the orbits of neighbouring Tethys and Dione from the point of view of dynamical stability. Large scale numerical simulations of 360 test particles within the coorbital regions of these four saturnian satellites were carried out for 4 × 10 5 yr or 1.6 × 10 8 revolutions of the innermost moon Mimas. The tidal forcing of the satellites' orbits was not taken into account in these simulations. We have quantitatively reproduced the Mimas-Tethys 4:2 and Enceladus-Dione 2:1 mean motion resonances in the system and devised a scheme by which the parameter space of the coorbital resonance is sampled uniformly by our test particles. We observe that 6 out of the 36 integrated horseshoe

particles of Enceladus escaped the coorbital region. All 54 tadpole particles remained stable. The main cause of instability for Enceladus coorbitals appears to be the overlap between the coorbital resonance and the 2:1 mean motion resonance between the particle and Dione. This leads particles with starting semimajor axes near the horseshoe-tadpole separatrix to be ejected from the resonance, as proposed by [START_REF] Morais | The effect of secular perturbations and mean motion resonances on Trojan dynamics[END_REF], over timescales of ∼ 8 × 10 7 revolutions of Enceladus. For Mimas we observe a larger number of coorbital escapes overall, both of tadpole (7/54) and horseshoe (29/36) librators. An analysis of the observed dynamical evolution

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 1 Introduction

The regular planetary satellite systems of the giant planets represent evolutionary end states of tidal evolution into or through orbital mean motion resonances. This occurred over time periods several orders of magnitude longer than the age of the solar system in terms of satellite orbital revolutions [START_REF] Peale | Orbital resonances in the solar system[END_REF][START_REF] Dermott | Dynamics of the uranian and saturnian satellite systems -a chaotic route to melting Miranda?[END_REF]. Recovering the initial state by inverting the dynamics in this case is not generally possible since the outcome of a resonant encounter under tides is probabilistic (Henrard and Lemaître, 1983). Nevertheless, an understanding of the dynamics of resonance capture has proven to be a powerful tool for constraining, not only likely initial configurations, but also the internal structure and geophysics of a planet and its satellites [START_REF] Peale | Orbital resonances, unusual configurations and exotic rotation states among planetary satellites[END_REF].

The trojan configuration, where a particle occupies the neighbourhood around either one of two stable equilibrium points in the Sun-planet-particle system leading or trailing the planet by 60 degrees in its orbit, are currently known to exist at Jupiter, Mars and Neptune [START_REF] Levison | The dispersion of the trojan asteroid swarm[END_REF][START_REF] Marzari | The MATROS project: Stability of Uranus and Neptune Trojans. The case of 2001 QR322[END_REF][START_REF] Scholl | Dynamics of Mars Trojans[END_REF]. These bodies appear to be stable in their current dynamical environment, and long-term dynamical studies point to lifespans of order the age of the solar system. Consequently, their existence at some planets and not others are clues to the prevailing conditions during solar system formation and early dynamical evolution

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT [START_REF] Morbidelli | Chaotic capture of Jupiter's Trojan asteroids in the early Solar System[END_REF].

It is then reasonable to assume that the same holds for satellite systems. Of all the planets, only Saturn is known to have satellites accompanied by trojans. The satellites Telesto and Calypso occupy the leading and trailing triangular equilibrium points of Tethys; the satellites Helene and the recently discovered Polydeuces [START_REF] Murray | S/2004 S 5: A new co-orbital companion of Dione[END_REF] play the same role at Dione. Although these bodies appear to be currently stable, the reasons why (a) only these two saturnian satellites are currently accompanied by trojans (b) there are no known trojans associated with the satellites of either Jupiter or Uranus, is still a matter of conjecture. Consequently, understanding what dynamical routes lead to long-term instability, as well as the timescale(s) over which these mechanisms operate, is necessary for explaining the current picture. In this paper we are concerned with what we call the "static" stability of saturnian satellite trojans, that is, stability over timescales within which tides do not evolve the orbits significantly. Although we expect tides to play an important role in answering the long-term stability question, it is known that sources of trojan instability exist in the conservative case [START_REF] Dermott | The dynamics of tadpole and horseshoe orbits. I -Theory[END_REF][START_REF] Morais | The effect of secular perturbations and mean motion resonances on Trojan dynamics[END_REF].

It is important to establish which of these operate in this system and to what extent. This will form a solid basis on which to study the effect of tides in a future paper.

In previous work, [START_REF] Burns | Dynamically depleted zones for Cassini's safe passage beyond Saturn's rings[END_REF] orbiting Saturn between the orbits of Mimas and Enceladus. In Figure S1 of their work, it appears that > 10 objects in horseshoe or tadpole orbits survive as such for this period of time. [START_REF] Mourão | On the stability of hypothetical satellites coorbital to Mimas and Enceladus[END_REF] specifically addressed the stability of particles in tadpole orbits with Mimas and Enceladus and concluded that, under the present satellite configuration, these are stable over 10 4 yr.

Here we are concerned with objects in either tadpole (ie trojans) or horseshoe orbits, which we refer to collectively as coorbitals. Furthermore, we limit our study to the Mimas-Enceladus-Tethys-Dione (or METD) system. The two inner satellites Mimas and Enceladus are dynamically linked to the outer, more massive Tethys and Dione respectively through mean motion resonances [START_REF] Greenberg | Orbital resonances among Saturn's satellites[END_REF]. We have carried out large scale numerical simulations of coorbital test particles perturbed by these satellites, over timescales more than an order of magnitude longer than previously attempted, to map out domains of long-term stability in phase space and their boundaries to regions where particles escape over the integration timespan.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Our results, while consistent with previous works, reveals that long-term instability does exist within the coorbital regions of Enceladus and, particularly, Mimas.

The paper is organised as follows: In Section 2, we discuss the dynamical model that we use as a background for this work. It is important to understand to what degree this set of satellites, or a subset of it, can be isolated from the whole, and still retain those dynamical features, namely the known mean motion resonances, which would enable a realistic study. In Section 3, we describe the setting up of the initial conditions of our test particles that we deem necessary for an even, unbiased, sampling of the coorbital phase space. In Section 4, we describe our results, concentrating on the dynamical sculpting of the coorbital zones of Mimas and Enceladus and the mechanisms responsible for the dynamical instability that we observe.

Satellite Model

Before one can proceed with large-scale numerical simulations of test particles, it is important to understand to what degree the model can represent, or reproduce, the actual dynamics. In addition, for reasons of computational expediency we aim to use the simplest model that can reproduce some main characteristics of the system. For this purpose we integrated several different satellite configurations: (a) six major satellites in order of radial distance from Saturn (ie Mimas-to-Titan or METDRT) along with the three "classical" trojans mentioned above, not including the recently-discovered Polydeuces (b) the Mimas-to-Dione four-satellite subset (METD) with the three trojans (c) the Mimas-Tethys system (MT) with the two Tethys trojans Telesto and Calypso (d) the Enceladus-Dione system (ED) together with the Dione trojan Helene. We employed a hybrid integration scheme available as part of the MERCURY package [START_REF] Chambers | A hybrid symplectic integrator that permits close encounters between massive bodies[END_REF]. This is an N-body symplectic code that switches to a Bulirsch-Stoer scheme during a close encounter between either a massive body and a particle or two massive bodies. The satellite masses and initial state vectors as well as the J 2 , J 4 and J 6 harmonic coefficients of Saturn's gravity field used were taken from [START_REF] Jacobson | The orbits of the major saturnian satellites and the gravity field of Saturn from spacecraft and Earth-based observations[END_REF]. We then carried a har-

The existence of numerous resonances among the saturnian satellites presents us

A C C E P T E D M A N U S C R I P T
ACCEPTED MANUSCRIPT monic decomposition of the numerical output using the Frequency Modified Fourier Transform (FMFT) algorithm developed by [START_REF] Šidlichovský | Frequency modified Fourier transform and its applications to asteroids[END_REF] implemented into an IDL program written by one of us [START_REF] Christou | The effect of near resonances on the secular precession of close orbits[END_REF].

The integration output was sampled with a period of 0.46875 days, slightly shorter than half the orbital period of Mimas to mitigate against frequency aliasing.

The results for the main satellite resonances are given in Table 1; those for the trojan librations are given in Table 2. The numerical estimates from the MET-DRT model give the best agreement with the observed amplitudes of libration of the Enceladus-Dione and Mimas-Tethys resonances of 1 • .5 and 97 • [START_REF] Allan | The evolution of the Mimas-Tethys commensurability[END_REF][START_REF] Sinclair | On the origin of commensurabilities amongst the satellites of Saturn[END_REF]. The same conclusion may be reached for the respective periods of libration [START_REF] Harper | The orbits of the major satellites of Saturn[END_REF] and forced actions (ie eccentricities and inclinations) corresponding to these two resonances. Note that the calculation of the libration amplitudes from Harper and Taylor's work is not a straightforward exercise. These authors provide the amplitude of the term containing the resonant frequency in the expression for the mean longitude of the satellites, not the libration amplitude of the critical angle itself.

The other models differ from METDRT in various degrees. METD differs by 10% and 1% respectively in the amplitude and period of libration of the φ e argument, while the differences for the φ I libration are 6% and 2%. The variation in the actions forced to the smaller satellites of the pairs by the resonances is found to be 1%. Most

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
of the resonant parameters for the two-satellite models MT and ED exhibit marked departures from these values. In particular, the amplitude of libration for those models is a factor of 7/2 larger for φ e and smaller by a factor of 2/3 for φ I in comparison to the other models.

Finally, we see that the parameters for the Trojan librations are relatively modelinsensitive in all cases. Our libration amplitudes for the Tethys trojans Telesto and Calypso are in good agreement with the values reported in [START_REF] Jacobson | Orbits of the saturnian satellites from Earthbased and Voyager observations[END_REF] and [START_REF] Porco | Cassini imaging science: initial results on Saturn's rings and small satellites[END_REF] of 1 • .2, 3 • .8 and 1 • .3, 3 • .6 respectively. They are somewhat smaller than those reported in Table 3.2 of Dermott and Murray (1999

) (2 • , 4 • )
and by [START_REF] Oberti | An upgraded theory for Helene, Telesto and Calypso[END_REF] (2 • .25, 5 • .25). The trojans' libration periods as reported by the latter authors are 6% shorter than ours. Our results for the Dione trojan Helene are in excellent agreement with the libration extrema reported by [START_REF] Greenberg | Orbital resonances among Saturn's satellites[END_REF] and [START_REF] Murray | S/2004 S 5: A new co-orbital companion of Dione[END_REF]. Our estimated libration period of 2.10 yr is bracketed by that reported by the latter authors (2.08 yr) and that in [START_REF] Oberti | An upgraded theory for Helene, Telesto and Calypso[END_REF] (2.14 yr).

One can draw several conclusion from these observations. Firstly, it appears that current knowledge of the saturnian system is sufficiently close to reality that simulating more complete configurations (ie more satellites, higher order planetary zonal harmonics) results in better agreement with the observed values of sensitive indicators of the system state. Indeed, during this process we became aware of a

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
scarcity of self-consistent numerical models of the resonant dynamics within multisatellite models such as METD or METDRT in the published literature; we expect the results reported here to serve as a useful reference in future work on the longterm dynamics of the system. More pertinent to this investigation is that the 4-satellite model (METD) appears to offer the sought-after compromise between computational expediency and model realism and is thus adopted as our base model.

Particle setup

Our system consists of the planet or primary body (mass M), a satellite or secondary body (mass m s ) with test particles in tadpole and horseshoe orbits (hereafter referred to as H and T respectively), and possibly other perturbing satellites. The particles' osculating elements are a (semimajor axis), λ (mean longitude), e (eccentricity), (longitude of pericentre), I (inclination) and Ω (longitude of ascending node).

Unsubscripted quantities refer to the test particles and quantities with subscript s refer to the secondary.

For the purposes of this work, we aim to place a parameter space grid of test particles within the coorbital regions of the satellites Mimas, Enceladus, Tethys and Dione.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

We use an initial condition setup based on current knowledge of the short-term and long-term stability of the three body problem. The short-term stability is reflected by our use of the semimajor axis cutoff (outermost horseshoe). The longterm is reflected by the use of classical secular perturbation theory and what it tells us about the forced eccentricity and inclination [START_REF] Morais | A secular theory for Trojan-type motion[END_REF][START_REF] Morais | Hamiltonian formulation of the secular theory for Trojantype motion[END_REF].

While unstable orbits in the three-body problem will remain unstable in the N-body problem, we test the stable orbits of the three-body problem under the perturbations of mean motion resonances (MMRs) as well as secular perturbations from additional satellites.

We choose to initially locate our test particles at the closest approach with the secondary as illustrated in Fig. 1. The expressions for the initial osculating semimajor axis and mean longitude for each particle are thus a = a s

(1)

λ = λ s + 60 • -l(i) ( 2 )
where l is a parameter that characterizes the amplitude of the coorbital which we choose as

l(i) = 1 - i N -1 60 • - 8(μ/3) 1/3 180 • 0.74 2 π (3) with μ = m s /M and i = 0, • • • , N -1. A parameterisation equivalent to Eqs. 1-3
uses the relative fractional semimajor axis a 0 = Δa/a s (Δa = aa s ) evaluated at [START_REF] Morais | A secular theory for Trojan-type motion[END_REF]. It is important to emphasize that a 0 and l are constants of the motion. The osculating quantities describing the particle's trajectory are Δa/a s and λλ s .

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT λ-λ s = 60 • (
The particles described by Eqs. 1-3 above are effectively spread at equal distances along an arc at the semimajor axis of the satellite. The factor 0.74 at i = 0 corresponds to the largest amplitude stable horseshoe, empirically derived from numerical experiments on the restricted circular three-body problem [START_REF] Dermott | The dynamics of tadpole and horseshoe orbits. I -Theory[END_REF]. This is appropriate here as the eccentricities and inclinations of the secondaries and test particles are not large compared to the mass parameter = (μ/3) 1/3 [START_REF] Namouni | Secular interactions of coorbiting objects[END_REF]. The choice of i = N -1 places a particle at the L 4 leading triangular equilibrium point. In this work we have used N = 10.

The corresponding expressions for the initial osculating eccentricity and inclination vectors e = e exp i and I = I exp iΩ are given by [START_REF] Morais | Hamiltonian formulation of the secular theory for Trojantype motion[END_REF] as e = e f + e p (4)

I = I f + I p (5)
where

e f = e s exp [i ( s + b(l))] ( 6
)
I f = I s exp [iΩ s ] ( 7 ) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
and e p = e p exp [i (γ(l)t + χ)] (8)

I p = I p exp [iΓ(l)t + Ξ] . (9)
The constants e p , χ, I p , Ξ are set by the initial conditions and the expressions for b(l), γ(l) and Γ(l) are given in [START_REF] Morais | A secular theory for Trojan-type motion[END_REF].

In the equations above, the "f" subscript denotes the forced component which depends on the orbital elements of the secondary and its position relative to the particle. The subscript "p" denotes the proper or free component which can be arbitrarily set to grid the initial conditions.

Fig. 2 shows a schematic representation of the eccentricity and inclination vectors defined in Eqs. 4-9, along with our chosen initial conditions. To obtain initial conditions for our simulations we decided to set the constants χ = b(l) and Ξ = 0.

Moreover, we set e p = 0, e s /2, e s and I p = 0, I s /2, I s so that e = (1 + j/2) e s exp [i ( s + b(l))] (10)

I = (1 + k/2) I s exp [iΩ s ] (11)
as illustrated in Fig. 2, with j, k = 0, 1, 2. Hereafter, individual particles will be referred to by using the values of i, j and k applicable to them, optionally preceded particles. We refer to population subsets (ie subset i = 3) using the i, j and k index notation and the starting values of the orbital elements in the main text except where we deem it preferable to use only one notation to avoid repetition.

It is pertinent to the ensuing discussions to note that the critical value l crit that marks the boundary between H and T motion (the H-T separatrix or HTS) is independent of μ. It is given by

l crit = π/3 -2 arcsin √ 2 -1 /2 36 • .094 (12)
and the corresponding a 0crit is equal to 8/3μ.

Results

Our simulations were run on a CONDOR mini-cluster comprising 4-6 high-spec PCs.

The integration timespan was 4 × 10 5 yr corresponding to 1.6 × 10 8 revolutions of Mimas or 5.6 × 10 7 revolutions of Dione. We consider that a particle escaped if its semimajor axis at the end of the integration is outside a band of half-width a s centered at that of the satellite and Δλ = λλ s is circulating. In the first instance,

we observe that all of the coorbitals of Tethys and Dione remained in libration. In contrast, 36 out of the 90 particles started as coorbitals of Mimas escaped, as did a small fraction of the Enceladus coorbitals (6 out of 90). In terms of libration type, all 6 Enceladus escapes were initially in horseshoe orbits, out of the 36 integrated

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
for that satellite. For the Mimas particles, 29 out of the 36 integrated horseshoes escaped, as did 7 out of the 54 tadpoles. In what follows we concentrate on the analysis of the cases of Mimas and Enceladus.

Enceladus

The 6 escaped Enceladus coorbitals all belong to the i = 3 bin, sharing an initial value of l (37 • .8) near that of the HTS (Fig. 3). [START_REF] Morais | The effect of secular perturbations and mean motion resonances on Trojan dynamics[END_REF]hereafter M00) found that the MMRs corresponding to the arguments 2λ Dioneλand 2λ Dione -λ-Dione (the e and e resonances according to the notation of Murray and Dermott (1999) which we adopt) are located deep into the coorbital zone at a 0 0 and can affect those coorbitals with a 0 a 0crit (or, equivalently, l l crit ) near the HTS. Those particles' coorbital period is sensitive to small perturbations particularly near the L3 equilibrium point at (Δa/a s = 0, λλ s = 180 • ). As a result, the time they spend within the MMRs varies from one coorbital libration to the next. As the e resonance also affects the eccentricity of the coorbital particle, the net result is a buildup in the eccentricity, which increases monotonically until the particle is ejected from the coorbital region by close encounters with either the secondary or nearby satellites. For the details of the mechanism we refer the reader to Fig. 6.4a and Section 6.4 of M00 (available at http://www.astro.mat.uc.pt/∼hmorais/thesis.ps.gz) but we stress that its efficiency depends not on only on the strength of a given

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
MMR but also its proximity to a 0 = 0.

The locations of these two MMRs at a 0 ∼ 0.05 and a 0 ∼ 0.17 for our model -which uses different starting conditions and parameters -agree closely with those of M00 (Fig. 4). It is also useful for completeness and the discussions to follow in the next section to note that the 4:2 MMR inclination triplet (I 2 , II and I 2 ) lies between l = 44 • and l = 55 • , near i = 0, 1, 2. As the inclinations of either Enceladus or Dione are small (∼ 0 • .02) the effect of these MMRs on the particles' stability should be negligible, consistent with our results.

The types of dynamical behaviour we observe (Fig. 5) are consistent with the proposed mechanism of M00: initially, the eccentricities increase, essentially monotonically, to values around 0.06 -0.07, typically within 25% of the integration time;

this is accompanied by moderate changes in the coorbital energy integral a 0 (Phase I) and also, something heretofore unreported in M00 and other dynamical analyses of this kind, erratic variation in I. Phase II is characterised by erratic variations in e and I due to moderately close approaches to Enceladus. Close approaches to nearby satellites require higher eccentricities and are not possible at this stage. Finally a close approach to Enceladus ejects the particle from the coorbital region (Phase III)

and the particle's eccentricity can now increase sufficiently to allow close encounters with the other satellites. Notably, the first two phases (ie erratic behaviour without escape) were also observed for particles adjacent to i = 3. This suggests that, given

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
enough time, those too will share the fate of their escaped siblings. The median escape time for all particles with i = 3 was ∼ 2.9 × 10 5 yr, 72% of the integration timespan, or 8 × 10 7 Enceladus revolutions.

Within the limitations imposed by small numbers statistics, there appears to be a preference for instability among particles starting with high proper eccentricity e p and low proper inclination I p . High-e p particles will encounter the secondary sooner than low-e p ones while high-I p provides a short-lived protection mechanism against close particle-secondary approaches.

Mimas

Mimas coorbitals, on the other hand, are more unstable. All escaping particles (36)

have turning points within 34 • of Mimas while tadpoles with values of l below that limit survive until the end of the integration. Overall, there is a trend towards faster escape as l approaches the secondary (Fig. 6). [START_REF] Burns | Dynamically depleted zones for Cassini's safe passage beyond Saturn's rings[END_REF] found a ratio of < 20/3 for H and T particles respectively after 10 3 yr. Our estimate for the same timescale is 31/54. The difference is not suprising since the two works sample particle initial conditions in different ways. escape times of several tens of percent of the integration timespan.

A closer, particle-by-particle inspection reveals that 9 out of the 13 particles belonging to the fast population escape within 1% of the start of the integrations or 4 × 10 3 yr; and that most of these particles (10 out of 13) belong to the i = 0, 1 set with l > 50 • . In addition, for the fast escapers only Phases II and III were observed;

the small, erratic variations in a 0 and monotonic increase in e, I that characterise Phase I are absent here, probably because they occur on timescales much larger than the observed timescale for escape. In contrast, the slow population is dominant for i = 3, 4 particles (or 32 • .0 < l < 38 • .4) while it is extant within the i = 2 set (l = 44 • .8) as well. Median escape time for the former set is 30 and 35 percent respectively of the integration timespan, 1.2-1.4×10 5 yr or 4.6-5.4×10 7 revolutions of Mimas. This is somewhat faster than for the corresponding slow escaping population of Enceladus. Within this population we see all three Phases that were found in the dynamical evolution of Enceladus coorbitals but not necessarily all for every particle.

Phase I-II transition typically takes place once an eccentricity of 0.03-0.04 has been reached. For particles with j = 1, 2 (44 • .8 < l < 51 • .2) this is close to the starting eccentricity so the particles can enter directly into Phase II at the beginning of the integration. We show such an example in Fig. 8 where particle 311 undergoes erratic, zero-average-slope variation in e and I characteristic of Phase II behaviour from the onset (top panels). We contrast this with the behaviour of particle 300 Apparent in the dynamical evolution of the two particles plotted in Fig. 9 is also a modulation of the semimajor axis, but not the longitude. This is observed in all coorbitals of Mimas perturbed by the 2:1 resonant potential of Tethys and is found to be of amplitude 45 -50 × 10 -5 in Δa/a s . We have re-examined the corresponding plots for the Enceladus particles and identified a similar modulation there, of amplitude between 10 and 25 ×10 -5 . The amplitude appears to be independent of l (with an exception at l l crit to be discussed shortly); it is probably due to the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
differential effect of the 2:1 MMR with Tethys on Mimas and the particle.

What is particularly relevant to this work is that the amplitude of this modulation exhibits a sudden increase at i = 3 (l = 38 • .6) ie in the vicinity of the HTS, reaching values (in units of 10 -5 ) of 70 for Enceladus and 62 for Mimas. From the phenomenological point of view, this indicates that a mechanism similar to that responsible for the observed instability of Enceladus coorbitals is also present in Mimas' case. The effect of a mean motion resonance, namely the 4:2 II MMR, on the semimajor axis is proportional to the time Δt that the particle spends within the MMR (Eq. 6.48 in M00) and this is maximized at i = 3 .

The larger extent of this mechanism at Mimas, causing instability for particles with values of i adjacent to i = 3, may be due to the I 2 and I 2 MMRs. In contrast to the Enceladus-Dione case, the inclination of Mimas and Tethys are not negligible (1 • .5 and 1 • .0 respectively) making these resonances comparatively strong.

In addition, they are located symmetrically with respect to II at a 0 = 0 (Fig. 9).

Hence, they are likely to affect large amplitude T or small-amplitude H motion for a good fraction of the coorbital libration period.

In parallel to what was observed at Enceladus, this mechanism affects the eccentricity as well as the inclination of the slow escapers. Regardless of the precise nature of the coupling between e and I, and whether it leads to actual growth in e (instead of zero-sum random jumps) and eventually escape from the coorbital zone,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

we have identified at least one escape mechanism for these Mimas coorbitals, related to the fast escapes for i = 0, 1.

These latter escapes are primarily caused by the 2:1 eccentricity-type MMR triplet (Fig. 9). In Fig. 10 we show an example of particle dynamics with i = 0 exhibiting slow circulation near the e resonance and corresponding slow variations in a and e. We posit that particles are transported in this way outside the region of stable H motion and escape; effectively these MMRs act to shrink the zone of coorbital stability around Mimas to just outside i = 2 (l > 44 • .8). The efficiency of this mechanism appears also to be dependent on the particles' starting eccentricity; the four surviving particles with i = 0, 1 (Fig. 6) have j < 2 (e < 0.03).

Escape for particles interior to i = 1 (l < 51 • .2) is facilitated by their initial random walk in a 0 caused by the interaction of motion near the HTS with the 4:2 inclination triplet as discussed earlier. In Fig. 11 we show an example of this mechanism in operation for a particle with i = 3.

The semimajor axis evolves until a sufficiently high value of a 0 is reachedtypically 0.001 -for the e resonance to start affecting the motion. The strength of this MMR is proportional to Tethys' eccentricity which is small (∼ 0.001), but not negligible compared to that of Dione or Enceladus. The particle then essentially becomes a fast escaper. Our confidence in this two-stage mechanism is bolstered by the observations that (a) all particles reach this value of a 0 just before they escape

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
(b) it explains naturally the gradual escape time profile in Fig. 6; deeper coorbitals need to random-walk their semimajor axis for a longer period before they achieve a high enough value to be "fed" into the region of fast escape.

Conclusions and Discussion

In this work we have investigated the stability of coorbitals of the saturnian satellites, in particular Mimas and Enceladus, in their current orbital configuration, over timescales for which the effect of tides is not significant.

All of the 54 Enceladus tadpole particles remained stable throughout the integration. A small fraction of our Enceladus horseshoe particles (6/36) escape through the mechanism proposed by [START_REF] Morais | The effect of secular perturbations and mean motion resonances on Trojan dynamics[END_REF]. This affects particles near the horseshoe-tadpole separatrix (HTS) where the interaction with a mean motion resonance (MMR) deep within the coorbital region results in eccentricity growth and escape. We thus confirm that the mechanism as studied by M00 in the Saturn-Enceladus-Dione-Particle problem is also efficient in a more complete satellite configuration with parameters closer to those of the actual Saturnian satellite system.

Surprisingly, these eccentricity-type resonances were also observed to affect the inclination of Enceladus coorbitals, a matter that deserves further investigation.

For Mimas, we have found a more complicated situation. 29 out of 36 horseshoe This qualitative difference between the two satellites with regards to the fast escaping particles we observed at Mimas has implications for the existence of small objects in orbit-sharing configurations with the two satellites. It implies, due to the time-reversibility of the conservative dynamics, that impact ejecta originating at Mimas can survive as large-amplitude horseshoe attendants of that moon for periods of 10 3 -10 4 yr with eccentricities no larger than 0.03 -0.04. If such objects exist, they may eventually be discovered by imaging surveys at Mimas' distance from Saturn.

Finally, small to moderate libration amplitude tadpoles of both satellites (l ≤ 36

• A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
at Enceladus and l ≤ 25 • at Mimas) survived until the end of the integration. If the absence of known librators is taken to represent the real picture, it implies that either such objects never existed or that they were ejected from the coorbital regions at some point in the satellites' history. The existence of tides in the Saturnian system offers fertile ground for investigating these questions further using tools such as the ones employed in this work. This will be the subject of a future paper. This is a fast escaper, leaving the coorbital zone after 6.5×10 3 yr. The two upper panels show the evolution of the semimajor axis and the eccentricty, while the bottom panel shows the critical angle of the e resonance ("ep" in Fig. 9). This angle exhibits circulation with a ∼ 10 yr period, and modulates both a and e. The longer period (∼ 50 yr)

apparent in e may be related to a similar period in a although there does seem to be some phase difference, particularly at the beginning of the integration. A yet slower modulation in e (several hundred years period) correlates well with a similar period in the inclination

  followed the dynamical evolution of an ensemble of particles orbiting Saturn among the major satellites over 10 3 yr. The A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT main motivation for their work was the risk to the Cassini spacecraft due to impacts with circumplanetary material. In the course of their study, they found that objects on tadpole or horseshoe orbits with the satellite Mimas and Titan in particular are stable over this timescale. More recently, Porco et al (2005) integrated the orbits of 1000 test particles for 4800 yr in the vicinity of Mimas' orbit in the context of their discovery of two small moons, Methone (S/2004S1) and Pallene (S/2004S2),

  measuring stick by which to gauge how "realistic" our model is. These resonances and their corresponding critical angles, are (i) the Enceladus-Dione eccentricity-type resonance with φ e = 2λ Dioneλ Enc -Enc (ii) the Mimas-Tethys inclination-type resonance with φ I = 4λ Tethys -2λ Mimas -Ω Tethys -Ω Mimas and (iii) the coorbital librations of Telesto and Calypso with Tethys (l r1 = λ Telλ Tethys , l r2 = λ Calλ Tethys ) and Helene with Dione (l r3 = λ Helλ Dione ). Here we use the standard notation for the orbital elements where λ is the mean longitude, Ω the longitude of the ascending node and the longitude of pericentre.

  by the name of the secondary eg ENCELADUS 302 or MIMAS 721. In the Figure captions we also provide the actual starting values of the orbital elements of the A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

  The escape time statistics are presented in Fig 7andshow two distinct populations of particles: "fast" escapers with typical escape times less than 10 4 yr, and a slowly-escaping population similar to that seen at Enceladus with characteristic

  Fig. 9 shows the locations of nearby multiples of the 2:1 MMRs with Tethys

  did 7 out of 54 tadpoles. Our results indicate the existence of two distinct, yet interacting, destabilising mechanisms operating on different timescales. Moderate-to-large amplitude horseshoes become quickly unstable owing to the positioning of Dione's 4:2 MMR eccentricity triplet just inside the radial boundary of Mimas' coorbital zone. The slower-acting mechanism, on the other hand, is reminiscent of the instability we observed for Enceladus coorbitals. It affects particles near the HTS and operates on a timescale comparable to the timespan of the integrations. In the absence of other plausible hypotheses, we attribute it to the symmetric positioning of Dione's 4:2 MMR inclination triplet within Mimas' coorbital region affecting near-HTS particles not just near L3 but also near the points of their maximum excursion in Δa/a s when λλ s ±60 • .

AA

  the test particles in (a, l) space where l = 60 • -(λ minλ s ). The particles are staggered at equal intervals on the a = a s circle from l = 0 (L 4 triangular equilibrium point) up to l = l max corresponding to the horseshoe stability boundary at a 0 = 0.74 derived empirically by[START_REF] Dermott | The dynamics of tadpole and horseshoe orbits. I -Theory[END_REF]. This parameterisation completely defines the energy integral of the coorbital resonance. 38 2 Top row: illustration of the decomposition of a coorbital particle's eccentricity and inclination vectors e exp i and I exp iΩ respectively into their free (or proper) and forced components. To first order, e forced e s and I forced I s . In the case of horseshoe libration, b(l) = 0. Bottom row: parameterisation of the initial e free and I free of the test particles in this work. . . . . . . . . . . . . . . . . . . . . . 39 Enceladus coorbitals as a function of angular distance from the L 4 triangular equilibrium point. The red curve denotes the number of escaped coorbitals per l bin (max. is 9 on the left vertical axis) while the blue curve is the median escape time per bin, replaced by 1.0 (right vertical axis) in bins where no particles have escaped. The top and bottom dashed horizontal blue lines mark times of 10 8 and 10 7 Enceladus orbital periods respectively. The dashed vertical line marks the tadpole-horseshoe separatrix. The shorter dotted lines denote, from left to right, the locations of e and e 2:1 MMRs with Dione. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4 The locations of the 6 different multiples of the 2:1 mean motion resonance with Dione within Enceladus' coorbital zone. These are indicated by horizontal segments on the left-hand side vertical axis. The two dashed horizontal lines indicate the outer boundary of the horseshoe region at 0.74 (upper line) and the maximum amplitude tadpole (lower line). Two test particles, one in a horseshoe (i = 0 or l = 56 • .7) and the other in a tadpole orbit (i = 4 or l = 31 • .5) have been superimposed to provide a sense of scale. . . . . . . . . . . . . . 41 As Fig. 3 but for Mimas coorbitals. The dashed blue vertical lines along the bottom indicate the positions of different mutiples of the 2:1 mean motion MMR with Tethys. From left to right: II , I 2 , I 2 , e , ee . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Escape statistics for Mimas and Enceladus coorbitals as a function of time. Bin size is 10 4 yr. . . . . . . . . . . . . . . . . . . . . . . . . . . Mimas coorbital particles over 3.1 × 10 5 yr. Particle 311 (l = 38 • .4, e = 0.03, I = 2 • .34; top) undergoes changes in a 0 causing it to alternate between tadpole and horseshoe libration while e, already high to start with, exhibits a flat, erratic variation characteristic of Phase II behaviour. Particle 300 on the other hand (e = 0.02, I = 1 • .56; bottom) starts with a smaller eccentricity which undergoes Phase I monotonic increase until values of 0.03-0.04 are reached. It leaves the coorbital zone at 1.8 × 10 5 yr (Phase III). Note how the evolution of I parallels that of e. . . . . . . . . . . . . . . . . 45 9 As Fig. 4 but for the 2:1 mean motion resonance with Tethys within Mimas' coorbital zone. Two test particles, one in a horseshoe (i = 0 or l = 57 • .6) and the other in a tadpole orbit (i = 4 or l = 32 • .0) have been superimposed to provide a sense of scale. . . . . . . . . . . 46 300 yr (or 10 5 Mimas revolutions), from the start of the integration, of particle MIMAS 312 (l = 38 • .4, e = 0.03, I = 3 • .12).

(A

  Figure 1: Christou, Namouni and Morais 2007, Saturnian satellite coorbitals

Figure 7 :

 7 Figure 7: Christou, Namouni and Morais 2007, Saturnian satellite coorbitals

Figure 10 :

 10 Figure 10: Christou, Namouni and Morais 2007, Saturnian satellite coorbitals

Table 1 :

 1 Dynamical parameters of the Mimas-Tethys and Enceladus-Dione resonances derived through FMFT decomposition of the following subsets of the Saturnian satellite system: (in order of increasing radial distance from Saturn) 6 satellites from Mimas to Titan (METDRT), 4 satellites from Mimas to Dione (METD) and the two satellite systems of Mimas-Tethys (MT) and Enceladus-Dione (ED). The four-satellite configuration was also integrated using a 15th order RADAU method (METDr) to test the robustness of the hybrid scheme used throughout the paper. The first and third lines provide the period (in years) and amplitude (in degrees) of the critical argument φ crit equal to φ e = 2λ Dλ E -E for the Enceladus-Dione resonance and φ I = 4λ T -2λ M -Ω T -Ω M for the Mimas-Tethys resonance. The second and fourth lines give the magnitude of the eccentricity/inclination forced on the smaller member of each pair (ie Mimas in MT and Enceladus in ED). Enceladus' forced eccentricity is given in units of 10 -5 while Mimas' forced inclination is given in degrees.

		METDRT	METD	METDr	MT	ED
	φ e	11.07 (1.505) 11.19 (1.663) 11.19 (1.661)		11.75 (5.37)
	e forced E	474	481	481		518
	φ I	69.59 (96.41) 68.30 (90.66) 68.20 (90.18) 63.89 (62.42)
	I forced M	1.573	1.572	1.572	1.564

Table 2 :

 2 Dynamical parameters of the three "classical" saturnian satellite trojans derived through FMFT decomposition of numerical simulations of the METD saturnian satellite model. In this Table, l r and a r are the libration maxima of λλ s and (aa s ) /a s respectively. The subscript 1 refers to Telesto (Tethys L 4 coorbital), 2 refers to Calypso (Tethys L 5 coorbital) and 3 refers to Helene (Dione L 4 coorbital). The first line for each object gives the libration period in yr, the second line gives the libration extrema of l r in degrees and the third line gives the libration extrema of a r in units of 10 -5 . System model notations are as in Table1.

		METDRT	METD	MT	ED
	l r1	1.905	1.905	1.905
		(58.7, 61.2)	(58.8, 61.1)	(59.1, 60.8)
	a r1	(-4.1, 4.3)	(-3.8, 4.0)	(-2.7, 2.8)
	l r2	1.907	1.907	1.907
		(-63.7, -56.4) (-64.0, -56.1) (-64.6, -55.6)
	a r2 (-11.9, 12.0)	(-12.6, 12.7)	(-14.1, 14.2)
	l r3	2.101	2.101		2.102
		(46.7, 76.6)	(46.7, 76.5)		(46.6, 76.7)
	a r3 (-62.0, 62.7)	(-61.8, 62.3)		(-62.5, 62.9)

Acknowledgements

Part of this work was carried out during a series of visits between the authors' respective institutions, funded by the following sources: Grant SFRH/BPD/19155/2004 (HM -Portugal), Observatoire de la Cote d'Azur (FN -France) and Grant PPA/V/S/1999/00628 (AC -UK).