A Noullez

K Tsiganis
email: tsiganis@auth.gr

Design of Low-Altitude Martian Orbits using Frequency Analysis

Keywords: Satellite Orbits, Mars Orbits, Frozen Orbits, Frequency Analysis

Nearly-circular Frozen Orbits (FOs) around axisymmetric bodies -or, quasi-circular Periodic Orbits (POs) around non-axisymmetric bodies -are of primary concern in the design of low-altitude survey missions. Here, we study very low-altitude orbits (down to 50 km) in a high-degree and order model of the Martian gravity field. We apply Prony's Frequency Analysis (FA) to characterize the time variation of their orbital elements by computing accurate quasi-periodic decompositions of the eccentricity and inclination vectors. An efficient, iterative filtering algorithm, previously applied to lunar orbiters, complements the method and is used to accurately compute the locations of POs/FOs, for a wide range of initial conditions. By defining the 'distance' of any orbit from the family of POs and using the relative amplitudes of the different components of the motion, we can build 'dynamical fate maps' that graphically depict the survivability of low-eccentricity, low-altitude orbits at every inclination, and can be used for efficient mission planning. While lowering the altitude generally enhances the effect of tesseral and sectorial gravity harmonics, we find this to have less consequence for low altitude Martian satellites, in contrast with the Lunar case. Hence, a high-degree ( 20 th ) axisymmetric model is adequate for preliminary mission design at moderate altitudes, but should be complemented at low altitudes by the methods described here. All families of POs and their spectral decompositions can be accurately and effectively computed by continuation in arbitrarily complex Martian gravity models, as our filtering algorithm requires only short integration arcs.

Introduction

Efficient design of satellite survey missions is typically based on a set of requirements and operational constraints. For example, the orbital geometry must be suitable for the mission at hand (e.g. low-altitude, near-circular, polar orbit for global surveys), while orbital variations due to the gravitational field of the target body (primary) be such that a reasonable lifetime is ensured, with minimal active control. The complexity of the gravitational field of realistic solar system objects does not necessarily facilitate this task.

Satellite motion around non-spherical bodies has been extensively studied, starting from the axisymmetric 'J 2 problem' (see [START_REF] Allan | The Critical Inclination problem : A simple treatment[END_REF]; [START_REF] Hughes | The 'Critical Inclination' : Another Look[END_REF]; [START_REF] Coffey | The critical inclination in artificial satellite theory[END_REF]; [START_REF] Jupp | The Critical Inclination problem -30 years of progress[END_REF]), where the only non-Keplerian potential term is the lowest-order zonal one. The general 'zonal problem' -in which the disturbing function contains only zonal harmonics of arbitrary order -can be cast into a one degree-of-freedom (d.o.f.) problem, if averaged over the mean motion of the satellite (i.e. ignoring possible tesseral resonances ; see [START_REF] Kaula | Theory of Satellite Geodesy[END_REF] and De Saedeleer (2005) for a closed-form theory), such that the mean semimajor axis a of the satellite becomes a constant of motion. The fixed points of this integrable dynamical system are the well-known Frozen Orbits (FOs), for which also the inclination i, the argument of pericenter (ω, or g in the usual notation of Delaunay variables) and the eccentricity e all become stationary. The right ascension of the ascending node (nodal longitude, Ω, or h in Delaunay form) precesses at constant rate. Irrespective of altitude, FOs have g = ±90 • and, for low eccentricities and high inclinations, e ∼ O(J 3 /2J 2 sin i), to a first approximation [START_REF] Coffey | Frozen orbits for satellites close to an Earth-like planet[END_REF]. For an extended discussion on higher degree/order approximations to FOs, the reader is referred to [START_REF] Gurfil | Motion near frozen orbits as a means for mitigating satellite relative drift[END_REF], [START_REF] Lara | Exploring Sensitivity of Orbital Dynamics with Respect to Model Truncation : The Frozen Orbits Approach[END_REF] or the textbook by [START_REF] Gurfil | Celestial Mechanics and Astrodynamics : Theory and Practice[END_REF]. For analytical models of FOs in Lunar orbiter problems, the reader is referred to the works of San-Juan et al. (2008) [START_REF] Abad | Analytical model to find frozen orbits for a lunar orbiter[END_REF] and Lara et al. (2009). [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF] also developed a canonical analytical theory for the high-order axisymmetric Lunar problem, while [START_REF] Delsate | Frozen orbits at high eccentricity and inclination : application to Mercury orbiter[END_REF] applied the theory to Mercurian satellites.

If the primary's potential has a relatively significant deviation from axial symmetry -as is the case for the Moon -it may be necessary to use models of higher complexity to adequately describe secular satellite motion. Non-axisymmetric potential terms add a second degree of freedom in the averaged problem and hence FOs with all four elements (a, e, i, g) being constant do not generally exist. Instead, Periodic Orbits (POs) of this 2-d.o.f. system can be found, for which i, e and g are not constant but suffer only low-amplitude, fast-periodic variations (see [START_REF] Tzirti | Effect of 3rd-degree gravity harmonics and Earth perturbations on lunar artificial satellite orbits[END_REF]), being however free from "slow", secular variations. Note that one can also define Frozen-onaverage orbits in tesseral potentials, removing h-periodic terms through a canonical transformation, as in [START_REF] Garfinkel | Tesseral Harmonic Perturbations of an Artificial Satellite[END_REF] for low-Earth orbits. As shown in [START_REF] Lara | Precise Analytical Computation of Frozen-Eccentricity, Low Earth Orbits in a Tesseral Potential[END_REF], inverting the generating function to re-introduce the shortperiod terms can be used to refine the nominal frozen solution, leading to an h-periodic orbit nearly identical (up to the order of the transformation) to the POs defined above.

Because of the near-preservation of orbital geometry and reduced excursions in orbital elements, periodic orbits are interesting for astrodynamics engineering and have thus received special attention (see e.g. [START_REF] Lara | On Periodic Polar Orbits in the Artificial Satellite Problem[END_REF] ; [START_REF] Lara | Searching for Repeating Ground Track Orbits : A Systematic Approach[END_REF] ; [START_REF] Lara | Repeat Ground Track Orbits of the Earth Tesseral Problem as Bifurcations of the Equatorial Family of Periodic Orbits[END_REF]; [START_REF] Tzirti | Effect of 3rd-degree gravity harmonics and Earth perturbations on lunar artificial satellite orbits[END_REF] and [START_REF] Lara | Design of long-lifetime lunar orbits : A hybrid approach[END_REF]). POs would appear as fixed points on a suitably chosen Poincaré section of the 4-d phase space of the averaged problem and they are in general non-isolated, i.e. they form 'families' that can be found by locating one such orbit and then smoothly varying some parameter (e.g. the inclination) of the system ; this is known as continuation method. A well-established method for computing families of period orbits is that of differential corrections [START_REF] Deprit | Natural families of periodic orbits[END_REF]). An efficient implementation algorithm can be found in [START_REF] Lara | On the numerical continuation of periodic orbits -an intrinsic, 3-dimensional, differential, predictor-corrector algorithm[END_REF]. The method also gives information on the linear stability of POs, and can be amended by differential algebra procedures (as in He at al. (2017)) for a more realistic propagation of uncertainties. This method has been successfully applied for finding repeat-groundtrack orbits in high-degree models of the Moon (Russel and [START_REF] Russell | Long-Lifetime Lunar Repeat Ground Track Orbits[END_REF] and the Earth (Lara and Russel, 2008).

Recently, in [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF] and [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF], we studied low-altitude satellite orbits around the Moon, using a novel frequency analysis (FA) algorithmthat we called Prony's method -as a tool for obtaining a global view of the secular dynamics, in arbitrary lunar gravity models. Our numerical FA method is quite efficient, as only short trajectory arcs are necessary for deriving an accurate quasi-periodic decomposition of the orbit; thus it can be applied to large sets of initial conditions. As found by different authors (Lara et al. (2009); [START_REF] Lara | Exploring the long-term dynamics of perturbed Keplerian motion in high degree potential fields[END_REF]) and confirmed by our own experiments [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF], a high-degree model is necessary to design low-altitude high-inclination orbits for the Moon, and a 'minimal' gravity model should be of degree at least l ≥ 7 for moderate altitudes and at least l ≥ 9 for low-altitude polar orbits. Moreover, going to lower and lower altitudes, the POs of non-axisymmetric models become increasingly displaced in (e, g) from the FOs of the corresponding axisymmetric problem, due to the relatively large C 2,2 and C 3,1 terms. The effect of these terms was essentially observed in [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF], who found that a satellite starting from a low-altitude FO (found in the axisymmetric model) can crash on the lunar surface, when subjected to a complete lunar gravity model, if no orbit maintenance is performed.

In [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF] we also presented an iterative filtering algorithm, as a tool that allows efficient computation of the families of POs of the 2-d.o.f problem. The algorithm is not really 'new'; it has been 'rediscovered' many times in celestial mechanics, for example by [START_REF] Couetdic | Dynamical stability analysis of the HD202206 system and constraints to the planetary orbits[END_REF], who used it to locate stable, resonant periodic orbits in two-planet systems, [START_REF] Noyelles | Expression of Cassini's third law for Callisto, and theory of its rotation[END_REF] and [START_REF] Robutel | Analytical description of physical librations of saturnian coorbital satellites Janus and Epimetheus[END_REF] who studied the 3:2 spin-orbit resonance problem and 1:1 resonant coorbital rotation, [START_REF] Dufey | Latitudinal librations of Mercury with a fluid core[END_REF] who analyzed the libration of Mercury, or Delsate (2011) who studied ground-track resonances around Vesta. All these authors used the NAFF method of [START_REF] Laskar | The chaotic behaviour of the solar system : A numerical estimate of the size of the chaotic zones[END_REF], for the frequency analysis part, which requires a trajectory arc that covers at least a few periods of the slowest-varying term. In [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF], we used Prony's FA method that only needs a short trajectory arc to fit all frequencies.

The main idea behind the filtering algorithm is to remove a selected frequency component from the decomposed eccentricity vector -in our case, the slow "proper" mode -to get to a pure Periodic Orbit PO that only contains 'fast' (forced) terms of low amplitude, overimposed on a quasi-stationary solution. These characteristics of the PO decomposition are explained by the theory of [START_REF] Wiesel | A Theory of Low Eccentricity Earth Satellite Motion[END_REF], which is devoted to the mathematical construction of near-circular periodic solutions and perturbations thereof, applied to reconstructing the trajectories of particular satellites. Despite the different approach and scope, Wiesel's elegant work also provides basis for the convergence of our algorithm : it demonstrates the existence of a single, slow 'proper' mode for near-periodic solutions and how, by subtracting this component from the solution (see Eq. (39) of that paper) one would be led ideally to a periodic solution (plus higher-order terms), with excursions that do not reduce to zero, in general.

Starting from an educated first guess (e.g. the FO of an axisymmetric model), we first compute the orbital evolution for a short time T ( 30 day ⊗ in the case of Mars, studied here, ⊗ denoting any parameter specific to Mars) to derive the quasi-periodic decomposition of the eccentricity vector. Then, we remove the slow (free precession) component, thus defining a new set of 'filtered' initial conditions. Repeating the integration for the new initial condition, we find a new orbit that is actually 'closer' to the PO, in the sense that the amplitude of the slow mode is now much smaller than before. As shown in [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF], typically after 2-3 iterations, the method converges to machine precision on the PO, for all initial eccentricities, inclinations and altitudes studied, such that the problem observed by [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF] goes away.

In the present paper we apply our FA method to analyze low-altitude orbits in axisymmetric and non-axisymmetric models of the Martian gravity field. Studies of 'special' (i.e. 2016) and references therein), so these orbits are still not as well-known as they are for the Earth or Moon and deserve study. We focus on using the iterative algorithm to identify the families of POs and FOs of each model and their differences. We are particularly interested in low-eccentricity POs (or FOs). The importance of such orbits for Mars can be understood by the fact that they are already employed in missions -e.g. the Mars Global Surveyor (i = 93 • , e ≈ 0.008, pericenter altitude 372 km). Analytical approaches to the study of their dynamics have been proposed by San-Juan et al. (2011), for an axisymmetric model. Here, we focus on the conditions needed to effectively retain very low altitude, namely a 'circular altitude' a -R ⊗ less than 100 km, where a is the semi-major axis of the satellite and R ⊗ the mean equatorial radius of Mars.

In the following sections we first describe some essential characteristics of our models of the Martian gravity field and recall the essentials of Prony's FA method and the filtering algorithm ; the reader is referred to [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF] and [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF] for more details. Then, we present our results on the computed families of POs/FOs and the dynamics in their neighborhood. Our conclusions are presented in the last section of this paper.

Martian gravity field and dynamical models

Let us examine whether the analytical and numerical tools that we developed for the study of orbits around the Moon (Tzirti et ) can be readily applied to the case of low-altitude Martian orbits. The two bodies in fact differ significantly, both in terms of their main physical properties, but also in the details of their gravitational field. That being said, the main ingredients needed to apply our formalism are similar, namely (a) wide-enough separation of the three principal frequencies of motion, which allow for smooth (i.e. non-chaotic) dynamics, and (b) availability of a very precise gravity field, which allows for an accurate propagation at low altitudes. For the altitude range considered here, between 0 and 500 km, the Martian gravity field dominates over other effects, like the solar three-body perturbation or the atmospheric pressure, which are neglected here.

Characteristic time scales

To study the secular motion of a low-altitude (massless) satellite around Mars, we need to average the disturbing function over the mean anomaly l. This is acceptable Averaging over the mean anomaly fast angle l of a satellite is allowed, provided the mean motion nrev(a) ≡ 2π/Trev(a) ≈ `GM a -3 ´1/2 (for non-spherical primaries, this Keplerian approximation to the mean motion should be corrected at least for the oblateness J 2 that slightly accelerates the revolution frequency ; see Murray and Dermott (1999) or Gurfil and Seidelmann (2016)) is much larger than the forced motion angular frequency 2π/Trot due to the primary's rotation. Also, to avoid resonances, the forced frequency and the slow free precession rate of the pericenter ġ = 2π/Ts(a) (computed here for e = i = 0 and accounting only for J 2 ) should be well separated. Both conditions are well satisfied for low-altitude orbits around all three primaries shown here. Eccentricity e F (J 2 &J 3 ; R, 90 • ) is the maximum eccentricity for the location of Frozen Orbits computed in the 'J 2 &J 3 problem' at i = 90 • and at the body surface R, while orbits with eccentricity larger than ec(a) at altitude a will collide with the body surface (note that ec(R) = 0). Also note the different signs of J 3 for the Moon and Mars with respect to the Earth.

when the mean motion, n rev (a) is fast enough and wellseparated from all other frequencies present in the problem : (i) the mean rotation rate 2π/T rot of the primary, which introduces a periodic forcing on the satellite orbit, and (ii) the (slow) precession rate of the pericenter argument g. Using the formulae of the averaged 'J 2 problem' (see e.g. Gurfil and Seidelmann (2016)) one gets ġ(a, e, i) = 3 J 2 n rev (a) 5 cos(2i

) + 3 8(1 -e 2 ) 2 R ⊗ a 2 (1)
which is as small as O(J 2 ). Table 1 lists the three characteristic periods for the Earth, Moon and Mars, and shows that the separation factor between the different frequencies is always in the range 12-14 for Mars. Note that the wide separation of the frequencies ensures that low-order resonances and chaotic phenomena are not present in the motion, and hence a meaningful quasi-periodic decomposition can be found by Frequency Analysis.

Martian gravity models and equations of motion

The Martian gravitational potential can be expressed as a spherical harmonics expansion in the body-fixed reference frame [START_REF] Kaula | Theory of Satellite Geodesy[END_REF])

V = - µ r ∞ l=0 R ⊗ r l × l m=0 P l,m (sin φ) [C l,m cos mλ + S l,m sin mλ] , (2)
where µ = GM ⊗ is the mass parameter and R ⊗ is the mean equatorial radius of Mars (3 397 km), (φ, λ) are the planetocentric latitude and east longitude respectively, r the distance of the satellite to the center of Mars, C l,m and S l,m the non-normalized adimensional gravitational coefficients, and P l,m the associated Legendre polynomials of degree l and order m. For m = 0 we get the zonal harmonic coefficients, for which we will use the notation J l = -C l,0 .

To integrate trajectories and use our technique to find lowaltitude orbits around a specific object, these coefficients have to be determined beforehand, either from the perturbations of the orbit of a relatively high-flying first satellite, or from a shape model coming from e.g. radar observations. In this work, we have made preliminary trials using different Mars gravity models of high degree and order, that is models MGS95J and MRO120D, coming from JPL (see [START_REF] Konopliv | A global solution for the Mars static and seasonal gravity, Mars orientation, Phobos and Deimos masses, and Mars ephemeris[END_REF], [START_REF] Konopliv | Mars high resolution gravity fields from MRO, Mars seasonal gravity, and other dynamical parameters[END_REF], [START_REF] Konopliv | An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data[END_REF] and references therein for a description), and Goddard Mars models 2B GMM-2B and model 3 GMM-3, described in [START_REF] Lemoine | An improved solution of the gravity field of Mars (GMM-2B) from Mars Global Surveyor[END_REF] and [START_REF] Genova | Seasonal and static gravity field of Mars from MGS, Mars Odyssey and MRO radio science[END_REF]. Of course, increasing (l, m) the model becomes more and more accurate, but also less and less tractable, as the number of terms increases quadratically with degree. Thus, as we did also in [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF], we tried to define a 'minimum model' for the purposes of this study, by progressively increasing the truncation degree l max from 8 to 24 and computing FOs and POs at a few selected inclinations with initial altitude a 0 = R ⊗ + 60 km. This is a procedure similar to the one used by [START_REF] Liu | Analytical investigations of quasicircular frozen orbits in the Martian gravity field[END_REF] for Mars orbits with a minimum altitude of 200 km. Systematic strategies to assess the minimal truncation order, needed to reach a given accuracy for different altitudes and inclinations were presented in [START_REF] Lara | Exploring Sensitivity of Orbital Dynamics with Respect to Model Truncation : The Frozen Orbits Approach[END_REF] for the Earth, and [START_REF] Lara | Exploring the long-term dynamics of perturbed Keplerian motion in high degree potential fields[END_REF] for the Moon, always based on progressively increasing the degree and checking the convergence of orbits to a given accuracy. We found the results to vary substantially up to l max ≈ 20 (much higher than l max = 9 used by [START_REF] Liu | Analytical investigations of quasicircular frozen orbits in the Martian gravity field[END_REF], as we are dealing with lower altitudes), beyond which the trajectories have very comparable secular evolutions ; also, for l max ≤ 20, all JPL and Goddard models gave similar results. Hence, we decided to adopt the Goddard Mars model 2B GMM-2B throughout this paper, truncated at l max = 20 and including all 0 ≤ m ≤ l terms, that we will now refer to as 20DM model. To quantify the effects of non-axisymmetric terms, we also employed a simpler 20ZM model, for which only the zonal terms J 2 -J 20 are retained to obtain an axisymmetric degree 20 model.

It is customary to employ the normalized gravity coefficients C l,m and S l,m , and fully normalized associated Legendre polynomials P l,m , to quantify the relative importance of terms in the above expansion [START_REF] Kaula | Theory of Satellite Geodesy[END_REF])

  P l,m C l,m S l,m   = (l -m)! (2l + 1) (2 -δ 0,m ) (l + m)! 1/2   P l,m C l,m S l,m   .
(3) Figure 1 is a graphic representation of the normalized gravitational coefficients J l,m ≡ (C l,m 2 + S l,m 2 ) 1/2 for both Moon and Mars. The relative strength of harmonics decreases rapidly with circular altitude a -R and degree l, as (R/a) l , following "Kaula's rule", i.e. J l,m ∼ l -2 . Such figures can be useful for comparing the relative importance of gravitational coefficients between two bodies, but also for assessing which (if any) can be neglected, to some approximation. For instance, for the Moon, J 5 is much smaller than the 'mean' and likely does not contribute as significantly to the dynamics as other coefficients of both higher and lower degree. Also, J 2 does not over-dominate the other terms, as C 2,2 is nearly as large.

Coefficients C 3,1 and J 7 appear on the other hand abnormally large and, as already noted by several authors (from [START_REF] Elipe | Frozen Orbits About the Moon[END_REF] to [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF]), they need to be taken into account for low-altitude lunar satellites. Note, however, that the figure should not be taken at face value, as different terms have different effects, depending on the geometry of the orbit considered and, in particular, higher-degree coefficients are more and more important at lower altitudes and larger inclinations. For example, the Moon's J 9 that falls right on the 'mean', controls the location of FOs at low-altitude polar orbits, as noted by (Lara et al., 2009) and [START_REF] Lara | Design of long-lifetime lunar orbits : A hybrid approach[END_REF] (this is also related to the very low values of J 11 , J 13 and J 15 , evident in figure 1). Gravitational coefficient amplitude, (

 C l,m 2 +  S l,m 2 ) 1/2 Mars, normalized adimensional gravitational coefficients  J l,m m = 0 m = 1 m = 2 m = 3 m = 4 m = 5 m = 6 m = 7 m = 8 m = 9
h P l m=0 J l,m 2 /(2l + 1) i 1/2
as a function of degree l, and decays approximately as l -2 (Kaula's rule, see [START_REF] Kaula | Theory of Satellite Geodesy[END_REF]).

For Mars, the situation is different, as J 2 dominates the dynamics, although not as much as for Earth. C 2,2 is relatively large -only 10 times smaller than J 2 -while it is 2 orders of magnitude smaller for Earth. In general, higher-degree terms and m = 0 terms are relatively stronger for Mars with respect to low-l and m = 0 terms, comparing to Earth. Hence, we expect the secular motion of low-altitude satellites around Mars, Earth and the Moon to be quite different.

In the Hamiltonian description, the motion of a satellite in the 20DM model takes the form

H n = H 0 + lmax l=2 H l,0 + lmax l=2 l-1 m=1 H l,m + lmax l=2 H l,l +H n⊗ , (4)
where H 0 = |u| 2 /2 -µ/r is the Keplerian part (u is the satellite's velocity), H l,0 , H l,m , H l,l describe respectively the zonal, tesseral and sectorial gravity harmonics and H n⊗ = -n rot H is added when working in a Mars-fixed rotating frame. The Hamiltonian is then expressed in canonical Delaunay variables (l, g, h, L, G, H), where

L = √ µa, G = L (1 -e 2
) and H = G cos i, and it is averaged over the mean anomaly l of the satellite in closed form. This procedure -which is valid for arbitrarily high eccentricity -is described in (De Saedeleer, 2005) and [START_REF] Tzirti | Effect of 3rd-degree gravity harmonics and Earth perturbations on lunar artificial satellite orbits[END_REF]. Other approaches for constructing the averaged Hamiltonian can be used and, in a recent work by [START_REF] Lara | Exploring the long-term dynamics of perturbed Keplerian motion in high degree potential fields[END_REF], the use of Kaulas' recursions [START_REF] Kaula | Theory of Satellite Geodesy[END_REF]) is pointed out as the most efficient method. Our symbolic code for this procedure is written in Mathematica and was checked against (De Saedeleer, 2006), [START_REF] Kaula | Theory of Satellite Geodesy[END_REF] and [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF]. The symbolic code produces FORTRAN expressions that are inserted into a propagation code, which gives the time evolution of (e, i, g, h) ; we re-mind that the semi-major axis a is constant in the averaged problem.

Frequency Analysis applied to Martian orbiters

Frequency Analysis is a powerful tool of time series analysis, and is particularly useful for planetary or satellite orbits, as it provides a quasi-periodic decomposition of the eccentricity and inclination vectors. A (real or complex) quasi-periodic "signal" u(t), is decomposed into a set of p periodic components

u(t) = p k=1 α k e i 2πν k t , (5)
with (real) frequencies ν k and (complex) amplitudes α k , using a sample of the signal over an "arc" of duration T , i.e. a finite number N of discrete values of u[n], separated by a constant time interval ∆t

u[n] ≡ u(n ∆t) n = 0, . . . , N -1 (6) = p k=1 α k ρ k n ρ k ≡ e i 2πν k ∆t . ( 7
)
Different ways of computing the spectrum of u[n] exist, as we described in detail in [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF]. The most common method is that of computing numerically the Fourier transform of u[n], assuming all frequencies ν k to be of the form ν k = k/T . The most refined version of this method is Laskar's NAFF algorithm [START_REF] Laskar | The chaotic behaviour of the solar system : A numerical estimate of the size of the chaotic zones[END_REF]) (see also [START_REF] Laskar | Frequency Map Analysis and quasiperiodic decompositions[END_REF] or [START_REF] Laskar | The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping[END_REF]). The main shortcoming of Fourier-based methods is that they are based on averaging the signal, thus requiring T to be (much) larger than the slowest mode that we need to resolve.

Prony's Frequency Analysis method

In [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF] we presented 'Prony's method', which is a fitting method ; this takes away the large-T limitation, so that short signal arcs are capable of capturing long-term oscillations. We repeat here only the essentials of the method and refer the reader to [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF] or [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF] for more details.

Prony's FA method is in fact the complex version of Prony's original method for real exponentials (de Prony, 1795) (see a modern description in [START_REF] Hamming | Numerical Methods for Scientists and Engineers[END_REF], (Hildebrand, 1987), [START_REF] Kay | Modern Spectral Estimation : Theory and Application[END_REF] or [START_REF] Noullez | Chaos Characterization in Hamiltonian Systems using Resonance Analysis[END_REF]), used to describe the time evolution of pressure in gas expansion. The main ingredient of the method is that any discretely sampled signal that is a simple sum of p exponentials as in (7) obeys a linear constant coefficients difference equation of order p u

[n] + a 1 u[n -1] + • • • + a p u[n -p] = 0 , (8)
and the p complex 'resonances' ρ k are the roots of the characteristic polynomial

ρ p + a 1 ρ p-1 + • • • + a p-1 ρ + a p = 0 (9)
of this difference equation. Prony realized that this gives a method to determine the frequencies, because finding the so-called prediction coefficients a j is a linear problem that we can solve exactly as soon as we have at least N ≥ 2p signal values. To do that, we write equation (8) for n = p, . . . , 2p -1 and solve the resulting p × p linear system

p j=1 u[n -j] a j = -u[n] n = p, . . . , 2p -1 (10)
for the p unknown a j , j = 1, . . . , p. Once the prediction coefficients are determined, ρ k , k = 1, . . . , p are found as the (complex) roots of the polynomial (9), using any standard algorithm (e.g. Newton's method in the complex plane), with a precision only limited by numerical accuracy (see [START_REF] Hamming | Numerical Methods for Scientists and Engineers[END_REF] or [START_REF] Hildebrand | Introduction to Numerical Analysis[END_REF]), and the frequencies ν k are then simply obtained as

ν k = Im {log (ρ k )} /2π ∆t = Arg (ρ k ) /2π ∆t (11)
because of the definition of ρ k eq. ( 7).

Once the frequencies ν k are known, the amplitudes α k , k = 1, . . . , p can be found easily, as the problem is linear in α k and we only need to solve another p × p linear system

p k=1 ρ k n α k = u[n] n = 0, . . . , p -1 . (12)
Even if the method is in theory exact -save for solving the polynomial (9) -the presence of noise would force us to use a least-squares Prony's method (see e.g. [START_REF] Kay | Modern Spectral Estimation : Theory and Application[END_REF]), and so N should be (much) larger than 2p, to allow accurate determination of p frequencies. Hence, the two linear problems described above become optimization problems that should be solved for increasing values of p, until the decomposition matches the signal to a pre-defined accuracy. Note that Prony's method assumes only distinct roots ρ k but, as we move in parameters space, distinct frequencies may collapse, or a frequency may become null (e.g. ġ at the critical inclination). Prony's method can be adapted to include secular drift terms and circumvent these issues, as described in [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF].

Spectral decomposition of orbits

We apply Prony's FA on the two complex variables u(t) ≡ e(t) × e i g(t)

(13)

v(t) ≡ sin i(t) × e i h(t) , (14)
i.e. the complex representations of the eccentricity and inclination vectors that determine secular motion. From the point of view of signal analysis, they have the advantage of being continuous when the g and h angles go through 2π ; also, the frequencies measured by FA are exactly related to ġ and ḣ, and their sign (prograde or retrograde) can be also determined, while analyzing e(t) or i(t) would prevent access to that information, as the spectrum of real variables is always symmetric (hermitian) between ν and -ν.

Orbital propagation is performed with an Adams predictor-corrector scheme, and the vectors are sampled at equidistant intervals ∆t; the sampling time has to satisfy the Nyquist criterion and thus must be chosen smaller than (half) the highest frequency. In the averaged problem, this is the rotational (forcing) frequency of Mars (i.e. the inverse of 1 Martian day), multiplied by l max . We employ a sampling time ∆t = 0.02 day ⊗ or 1/50 of a Martian day 1 .

The total length of the signal should be N 2p, and also such that the slowest mode must have had the time to vary 'significantly'. In practice, as we found in [START_REF] Tzirti | Secular dynamics of a lunar orbiter : a global exploration using Prony's frequency analysis[END_REF], we need N to be at least about 1/8 th of the slowest period. We decided to employ T = 30 day ⊗ (i.e. N = 1500 points), which allows reaching machine precision for all frequencies and amplitudes, down to ν ≈ 0.004 day ⊗ -1 for nearly all initial values of (i 0 , e 0 ) at given initial circular altitude a 0 -R ⊗ . Still, our trajectory arc and data length N is at least an order of magnitude shorter than what traditional Fourier-based methods would need, for the same accuracy.

Global view of low-altitude dynamics

We wish to obtain a global view of the dynamics, building a frequencies and amplitudes 'map', for a wide range of initial conditions, at sufficient resolution. We restrict ourselves to low altitudes, down to a few tens of kilometers. Performing test integrations, we found that the altitude range 50-500 km has similar dynamics, apart from the percentage of colliding orbits at low altitudes a 0 -R ⊗ . Hence, we present here the complete results only for a 0 = R ⊗ + 60 km. All orbits have h 0 = 0; this only controls the initial phase of the fast periodic component on the orbit. In a similar way, g 0 controls the initial phase of u. As we want this to be close to a FO/PO, we set it equal to g 0 = -90 • . As shown in [START_REF] Delsate | Frozen orbits at high eccentricity and inclination : application to Mercury orbiter[END_REF], the position of FOs is controlled by the odd-order zonal harmonics, and J 3 for Mars has the same sign (positive) as for the Moon, opposite to Earth.

Our FA maps are obtained for orbits starting on a i 0 -e 0 grid, spanning the range 0 < i 0 ≤ 90 • and 0 < e 0 < e c , where

e c ≡ a 0 -R ⊗ a 0 ( 15
)
is the collision eccentricity for circular altitude a 0 -R ⊗ , for which the satellite will crash on Mars surface. The quasi-periodic decompositions of u(t) and v(t) take the form (in order of increasing frequencies)

u(t) = α 0 + α s e i 2πνs t + k =0 α k e i 2πkν1 t (16) v(t) = k =0 β k e i 2πkν1 t ≈ sin(i 0 ) e i 2πν1 t (17)
with basic frequencies ν 0 ≡ 0 (constant), ν s and ν 1 . An example of the decomposition of the eccentricity vector u(t) is given in figure 2, with the first constant term α 0 shown as the thick black star ; the second term of ( 16) is displayed as the dotted blue curve, describing a circle of radius α s in the complex plane, traveled along at constant angular velocity 2πν s . The sum in the third term of ( 16) is represented by the thick red curve that closes itself after a period 1/ν 1 . The full eccentricity vector u(t) is the sum of all three terms and is shown as the thin black line. with full black circles separated in time by one Martian day. The two frequencies ν s and ν 1 are in general incommensurate, so the curve for u(t) doesn't close itself. Terms related to the 'fast' frequency ν 1 are present in both u(t) and v(t). These are forced terms, related to the rotation of Mars and the (slow) regression of h (caused mainly by the J 2 term at small inclinations), whose rate in the 'J 2 problem' is

ν 1 = 1 2π ḣ(a, e, i) -n rot (18) ḣ(a, e, i) ≈ -3 2 J 2 n rev (a) cos(i) (1 -e 2 ) 2 R ⊗ a 2 .(19) Note that ν 1 ≈ -1 day ⊗ -1
, accelerates only by ≈ 4% as i 0 goes from i 0 = 90 • to i 0 = 0 (see figure 3) and depends only mildly on e 0 . Hence, it is practically a function of i 0 only, for low-altitude Martian orbits. For v(t), only the single frequency ν 1 is present, the amplitude stays nearly -0.0005 0 0.0005 e × cos(g) -0.0095 -0.0090 -0.0085 e × sin(g) a = R ⊗ + 60 km, g 0 = -90°, h 0 = 0, i 0 = 90°, e 0 = 0.00907, t = 10:20 Figure 2: Graphic representation of the frequency decomposition of the eccentricity vector u = e cos(g) + i e sin(g) for a Martian polar orbit with i 0 = 90 • , e 0 = 0.009, g 0 = -90 • and h 0 = 0, at a 0 = R ⊗ + 60 km. The thin black line is the 30-days integration, and the full black circles, attached to this line, are separated by a time interval of 1 day ⊗ , showing stroboscopically the slow eccentricity drift. The black star denotes the constant found in the spectral decomposition, corresponding to the Frozen Orbit α 0 ≡ i e F sin g F , the dotted blue circles corresponds to the (extrapolated with empty circles) proper motion of g (the single 'slow' circular term of frequency νs), while the thick red line shows the sum of all the forced terms of the decomposition of frequencies kν 1 . The sum of these three components gives a very accurate reconstruction of the black line. Figure 3: Forced motion frequency ν 1 of the g and h angles for the 20DM (squares) and 20ZM model (diamonds) for g 0 = -90 • , h 0 = 0, e 0 = 0.00252 and an altitude 60,km, as a function of the inclination i 0 . The periodic frequency of both angles is the and is independent of g 0 and h 0 , and it is exactly equal to -nrot/2π ≡ -1.0 day ⊗ -1 for i 0 = 90 • . The dashed curve is the frequency -nrot/2π corrected by the effect of the J 2 coefficient (eq. ( 19)), that gives the main contribution to the slight retrograde motion of the node longitude h. Note also that for Mars, the tesseral and sectoral coefficients have no apparent effect on the frequency ν 1 . Figure 4: Proper motion frequency νs for the g angle for the 20DM (circles) and 20ZM model (diamonds) for g 0 = -90 • , h 0 = 0, e 0 = 0.00252 and an altitude 60 km, as a function of the inclination i 0 . This slow frequency is almost independent of g 0 and h 0 , has a period that is at least 12 day ⊗ (νs < 0.085 day ⊗ -1 ) at that altitude, and goes to zero around i 0 ∼ ic ≈ 63 • . The dashed line is the frequency that would be due to the effect of the J 2 coefficient only at that altitude (eq. ( 1)), and one can see that it gives the main contribution to the proper frequency νs. Once again, the axisymmetric and full model frequencies are undistinguishable.

constant at its initial value sin(i 0 ) and h circulates freely at constant rate. The fast dynamics of (h, i) is thus essentially unaffected by the motion of the pericenter (g, e), but it mildly affects the latter through the non-axisymmetric terms of the gravity field of the rotating primary.

The most significant eccentricity variations are due to the second term in (16) that corresponds to the slow proper frequency ν s that, as we can see in figure 4 is indeed well separated from n rot /2π 15 and from ν 1 1, taking values in the range [-0.020, 0.085] day ⊗ -1 as a function of i 0 , and is very close to ġ/2π as computed by Eq.( 1). The slight deviations seen at small inclinations are primarily caused by the C 2,2 term. Again, ν s depends only mildly on e 0 . However, the corresponding amplitude α s depends strongly on both i 0 and e 0 , as shown in figure 5.

The constant term α 0 of frequency ν 0 ≡ 0 corresponds to the mean position of u(t), around which the orbit is librating. Its value is, in general, not equal to the mean of the signal u over T , since the other terms would average to zero only if T is an integer multiple of both 1/ν s ≡ T s and 1/ν 1 ≈ T rot . We also stress that it is only thanks to the fitting properties of Prony's FA that this term can be determined accurately, using only a short orbit arc. For all orbits where this term is non-negligible, the amplitude α 0 is found purely imaginary, with complex phase Arg(α 0 ) = ±90 • . Hence, if we write α 0 ≡ i e F sin g F , the mean position of the pericenter is g F = ±90 • , as for Frozen Orbits in the axisymmetric problem. In fact, (e F , g F ) are almost exactly those of the FOs of 20ZM, as we verified -both models give the same α 0 and the orbits freeze in eccentricity e and pericenter g, when integrated in 20ZM. Of course, Figure 5: Slow mode amplitude |αs| as a 2-d grey levels map in the plane of initial conditions i 0 -e 0 , with fixed initial orientation angles g 0 = -90 • and h 0 = 0. White zones indicate orbits colliding in less than 30 days. Superposed on top are the location of Frozen Orbits FOs eccentricity e F ≡ |α 0 | found by FA for which the pericenter argument g F is -90 • (white circles) or +90 • (yellow stars). The regions with minimum slow mode amplitude |αs| closely follows the FO with the same orientation, and the slow mode amplitude increases linearly when moving away from the FO with the same direction, and even quicker when it has the opposite direction. This leads to collisions with the Martian surface at high eccentricity for i 0 < 13 • , all eccentricities around the critical angle i 0 ∼ 63 • , and at low eccentricities for i 0 > 75 • . these initial conditions do not 'freeze' exactly in 20DM (as we will see in section 5). Again, α 0 has strong variations with i 0 (circles and stars in figure 5, and squares in the left of figure 6), but is found to depend only mildly on e 0 (right of figure 6).

In figure 5, one can see that, when g F and g 0 have the same sign, the amplitude |α s | isolines closely follow the variations in position of the FO e F sin g F , which denotes both the mean position of the 2-d.o.f. PO and the mean position of the full orbit, as both slow and fast modes are periodic with zero mean. The slow motion amplitude is found to increase linearly like δe F ≡ |e 0 -e F | at any given inclination i 0 . In fact, the proper amplitude |α s | can be considered as the 'distance' of an orbit from the PO at the same inclination. Indeed, for vanishing |α s |, the motion becomes strictly periodic at frequency ν 1 for POs ('frozen' for FOs in 20ZM). The fact that |α s | ∼ δe F is probably a consequence of the large separation between the frequencies ν s and ν 1 that leads to nearly decoupled evolution, as noted in the previous section. Still, u(t) has a non-zero periodic component that slightly separates the true POs from the FOs (see figure 7), and complicates the search for initial conditions for POs (section 5).

Figure 6 shows a graph of e F ≡ |α 0 |, |α s | and α p for given initial eccentricity e 0 (shown as the thin horizontal line on the left) as a function of initial inclination i 0 (left), and for given initial inclination i 0 = 90 • , as a function of initial eccentricity e 0 vertically (right), both with g 0 = -90 • . The amplitudes |α k | of all modes ±k ν 1 for u(t) are found to vary only slightly with e 0 (see right of figure 6), depending mostly on i 0 , but the total forced amplitude α p defined as α 2 p ≡ k =0 |α k | 2 does not vary much with i 0 . Though α p is generally small compared to the other two terms, it never goes to zero, so that FOs with constant orbital elements do not exist in a non-axisymmetric gravity model and POs are frozen in the strict sense only when averaged over the period 1/ν 1 (see also [START_REF] Lara | Precise Analytical Computation of Frozen-Eccentricity, Low Earth Orbits in a Tesseral Potential[END_REF]). Note that ν 1 and ν s are both practically identical between the full (20DM) and the corresponding 20ZM axisymmetric model.

In many applications, the amplitudes |α 0 |, |α s | and α p are more interesting than the frequencies : they correspond to the mean eccentricity (for |α 0 |), the proper eccentricity 2 |α s | and a 'fast' periodic kick α p . Depending on the relative sizes of these amplitudes, the orbit will be circulating in g if |α 0 | < |α s | + α p , or librating otherwise. If g 0 has the same sign as g F , |α s | + |α 0 | ≈ e 0 , so that |α s | ≈ δe F . This implies that the amplitude |α s | decreases when going towards and vanishes at the PO, 'pointing' to its location (see figure 6 (right)). Thus, Prony's FA method provides a measure of orbital 'distance' from the relevant PO.

The magnitude |u(t)| is equal to the instantaneous eccentricity e(t). If e(t c ) ≈ e c for some time t c , the satellite crashes on the Martian surface ; this occurs when the sum in ( 16) equals e c but, because the frequencies ν s and ν 1 are well separated and non-commensurate, it suffices that e max ≡ |α 0 | + |α s | + α p ≈ e c for the satellite to crash in a time at most T s /2 + T rot (half a secular cycle). 2 As it is usually called in celestial mechanics.

Given the above discussion, it is clear that the orbit can crash either if e F or δe F are large, which can be both true for e 0 ≈ 0. This is shown graphically for polar orbits with a 0 = R ⊗ + 60 km in the right of figure 6. The orbits crash both for e 0 > 0.0163 (less than the collision eccentricity e c = 0.0174) or e 0 < 0.0011 (see also figure 5 on how low-e orbits crash for 75 • < i 0 ≤ 90 • ). This is what [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF] had noticed for low lunar orbiters and is the reason why a high-degree gravity model and knowing the location of the FOs are necessary for these studies.

Locating FOs/POs by filtering

As noted in the previous section, the periodic amplitude α p never vanishes for a non-axisymmetric gravity model, so Frozen Orbits with strictly fixed elements, free from periodic oscillations at frequency ν 1 , do not exist. However, Periodic Orbits have proper amplitude |α s | = 0 and will only have very small eccentricity variations (i.e. of size ∼ α p 10 -3 ) and so will have minimum eccentricity and altitude variations for a given i 0 . If we are moreover interested in minimizing altitude variations whatever i 0 , we can search for POs with e F ≈ 0. Such POs would have only fast periodic oscillations (amplitudes α k ) around a zero mean eccentricity ; the pericenter argument will certainly circulate, but the orbit will be quasi-circular.

In this section, we apply the iterative filtering algorithm presented in [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF] to compute the -0.0005 0 0.0005 e × cos(g) -0.0095 -0.0090 e × sin(g) a = R ⊗ + 60 km, g 0 = -90°, h 0 = 0, i 0 = 90°, e 0 = 0.00873, 20ZM/20DM Figure 7: FA decomposition of orbits (thin black lines) started exactly at the same initial conditions for the axisymmetric 20ZM (diamonds) and the full 20DM model (circles). Both orbits share exactly the same FO location (shown as a black star), but the axisymmetric orbit only contains a slow proper motion (indigo diamonds) that precisely describes the orbit itself. The full-model orbit also contains a periodic component (red curve obtained by FA) whose instantaneous value has to be subtracted from the initial conditions, to derive the proper motion (blue circles). The PO is found almost exactly if we remove the slow component from the initial orbit, to derive a periodic part that nearly coincides with the PO and whose initial condition is the red square).

locations and spectral decompositions of POs for all values i 0 at given a 0 -R ⊗ . We recall that locating a PO for a given i 0 implies finding (e 0 , g 0 ) points in a suitable h 0 section (we took h 0 = 0) where the real and imaginary parts of α s vanish simultaneously. As the complex amplitudes are smooth functions of the initial conditions and we have two 'free' parameters (e 0 and g 0 ) available to cancel the two components of the amplitude α s , one can expect to find a single periodic solution, for each i 0 and a 0 -R ⊗ . Note that ν s stays constant as α s goes to 0 when changing e 0 and g 0 , unlike in the 'J 2 problem' where the orbit freezes not because the proper amplitude goes to zero, but because the frequency ν s vanishes at the critical inclination i c , so that FOs exist for every e 0 at that particular inclination i c -however, g must be an integer multiple of 90 • , as becomes evident when considering the problem at order O(J 2 2 ) or higher ; see e.g. [START_REF] Coffey | The critical inclination in artificial satellite theory[END_REF], [START_REF] Broucke | Numerical integration of periodic orbits in the main problem of artificial satellite theory[END_REF], [START_REF] Gurfil | Motion near frozen orbits as a means for mitigating satellite relative drift[END_REF] and [START_REF] Gurfil | Celestial Mechanics and Astrodynamics : Theory and Practice[END_REF].

Iterative filtering algorithm

One might naively think that the mean orbit position α 0 ≡ i e F sin g F found by Prony's FA gives the initial conditions (e F , g F ) of the PO but this is not the case, as figure 7 demonstrates. In fact, α 0 is the position of the FO that can indeed be computed immediately by Prony's FA, in axisymmetric models. In 2-d.o.f. models, it is necessary to separate the slow and periodic components. Thus, the PO is found only after 'filtering' (i.e. removing) the slow component. But the periodic part has non-zero amplitude on the PO, and is actually needed for computing the PO's initial condition. This is the basis of the 'iterative filtering algorithm' that we will now describe briefly (see [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF] for details).

As described also in the Introduction, we start the procedure by adopting a first guess u (0) (0), v (0) (0) = {e 0 exp (i g 0 ) , sin(i 0 ) exp (i h 0 )}, for which we compute its time evolution for a short integration time T (30 day ⊗ ) and derive its quasi-periodic decompositions ( 16) and (17), using Prony's FA method. Then, we remove the slow component α s from the reconstructed time series, thus obtaining a 'filtered' orbit

{ u (t), v (t)} ≡ u(t) -α s e i 2πνs t , v(t) (20)
=    k α k e i 2πkν1 t , k =0 β k e i 2πkν1 t    ,
which is by construction periodic with frequency ν 1 . This operation does not change v(t), i.e. v (t) = v(t). Hence, the 'new orbit' has the same (i 0 , h 0 ) as the original.

If the α 0 and α k 's were completely independent of (g, e), the filtered orbit { u (t), v (t)} would be precisely the PO. This is however not true, as the filtered orbit is generally not an exact solution of the equations of motion, for these initial conditions { u (0), v (0)}. But, given that the main frequencies ν s and ν 1 are well-separated and no linear combinations with significant amplitude appear in u(t) and v(t), the proper mode defines a circle on the (e cos g, e sin g) plane, centered around α 0 . Elimination of the slow term sets the 'new' orbit at the center of this circle, adding the small periodic component. Thus, the 'new' orbit will have a slow part that is smaller than the original one and closer to the PO.

Using the new initial conditions, we iterate the procedure of integrating and filtering. Here, as in [START_REF] Noullez | Satellite orbits design using frequency analysis[END_REF] for the Moon, the procedure always converged to a PO in a matter of r = 2 -3 iterations, with a nearexponential decrease of |α s |. However, there is no rigorous proof that this algorithm always converges, despite attempts e.g. by Noyelles et al. (2012). Nevertheless, we did not encounter any case where the algorithm failed. Note also that the computational time for applying Prony's FA and filtering is negligible with respect to the integration time of the orbit.

Low-altitude Martian POs

Using our iterative filtering algorithm, we computed the POs (g p , e p ) for an initial altitude a 0 = R ⊗ + 60 km, and for 0 < i 0 ≤ 120 • . As starting guess we used the FO (e F , g F ) found by Prony's FA. The algorithm converged after at most 2-3 iterations, with a residual amplitude |α s | ≤ 10 -8 . The only exceptions to this fast convergence were found near the critical inclination i c ≈ 63.4 • , where the POs computed eccentricity becomes larger than e c so that all these orbits would crash on Mars surface. As in (Lara

a = R ⊗ + 60 km, h 0 = 0, 20DM(3) ℑ(α 0 ) ± α p FO ℑ(α 0 ) PO e 0 sin(g 0 ) PO ± α p FO J 2 +J 3
Figure 8: Initial (signed) eccentricity ep sin gp of the POs found after 3 iterations of the filtering algorithm (blue circles), for 0 < i 0 ≤ 120 • (note the symmetry around polar orbits) and a 0 = R ⊗ + 60 km. The modulus of the constant term ±|α 0 | that gives the mean position e F sin g F (i.e. FOs of 20ZM) is also shown as a black dashed line. The turquoise band shows the amplitude ±αp of the fast, forced motion (also shown as red diamonds around the e = 0 axis). The true POs are always close to the FOs (inside the periodic band), but their initial conditions differ due to the periodic part. The dashed magenta line is the locus of FOs in the simple 'J 2 &J 3 problem', which apparently cannot serve better than a good starting point for the algorithm, as it can be far from the true POs (see also figures 9 and 10). et al., 2009), we display the families of POs using e F sin g F instead of e F , to avoid an apparent slope discontinuity when g F jumps between ±90 • and ∓90 • as e F goes through zero. For POs, the initial pericenter argument g p is not exactly ±90 • and depends slightly on the initial nodal longitude h 0 . However, the amplitude of the periodic motion is small, so that g p will stay close to ±90 • in most cases.

Figure 8 shows that indeed e p sin g p and e F sin g F are smooth functions of i 0 and POs are always close to the corresponding FOs, but are not the same. In fact, the PO has variations of width ±α p around the FO. The graph also shows that the eccentricity of low-altitude Martian POs will always be small. The discontinuities seen around the critical inclinations i c ≈ 63 • or 117 • are real, and correspond to inclinations where the PO's eccentricity grows dramatically and goes above e c , so that non-collisional orbits cease to exist in these regions. POs with very small mean eccentricity exist for i 0 ≈ 13 • and i 0 ≈ 66 • or 114 • at 60 km altitude. As we will see in section 6.2, these inclinations allow to reach even lower altitudes. The PO's pericenter argument is always close to ±90 • , except for very low inclinations, where g 0 is determined only by the periodic component and is ill-defined, as e → 0. POs generally have g = -90 • , except for 0 < i 0 < 13 • and around i 0 = 66 • or i 0 = 114 • ; these are the only inclinations at this altitude for which we can get stable POs with g in the Northern hemisphere.

a = R ⊗ + 60 km, g 0 = -90°, h 0 = 0, i 0 = 90° | 60°, T = 23.5
Figure 9: Orbital evolution (thin black lines) for i 0 = 90 • (circles) and i 0 = 60 • (diamonds), started at e 0 = 0.0008 and e 0 = 0.00001 respectively, and a 0 = R ⊗ + 60 km. Both initially quasi-circular orbits eventually reach ec (violet, dashed circle) and crash on the Martian surface after 23.5 day ⊗ and 63 day ⊗ respectively, due to the slow rotation of the eccentricity vector (blue and indigo symbols) around the respective FOs (star and cross) and POs (red and brown curves). For the polar orbit, it is the small periodic component that actually causes the crash. The crash for i 0 = 60 • is detected using an integration up to T = 30 day ⊗ , thanks to the fitting/extrapolation properties of Prony's method. An orbit started at the location of the FO of the 'J 2 &J 3 problem' (magenta curve) at i 0 = 60 • only barely avoids crashing.

Dynamics near low-altitude POs

POs (and FOs) are interesting for practical purposes, as they dictate the dynamics in their phase-space neighborhood. For low altitudes in particular, the location of a PO can cause a quasi-circular orbit to crash on the primary's surface. This unobvious effect was already noticed by [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF], who studied the crash of an Apollo 16 subsatellite on the Moon after 34 days, while the nominal mission duration for 100 km altitude was one year. The crash was caused by a continuous growth of the subsatellite's eccentricity, that [START_REF] Knežević | Orbit maintenance of a lunar polar orbiter[END_REF] showed originating from l > 3-degree gravity terms that result into a slow rotation of the eccentricity vector around the corresponding FO. Hence, as shown below, having accurate positions of FOs or POs, we can build 'dynamical fate' maps for all initial conditions in their neighborhood.

Dynamics around FOs/POs

Orbits started at (e p , g p ) (or (e F , g F ) for 20ZM) have zero proper amplitude α s while, for any other initial condition e 0 , g 0 ), they have e 0 > e p , we will then have e(t) slowly oscillating within e 0 and 2e p -e 0 < e 0 , with g circulating if e 0 -e p > e p and librating otherwise. On the other hand, for e 0 < e p , the slow rotation around the PO lifts the eccentrity up to 2e p -e 0 > e 0 ; the satellite will crash if this is larger than e c . Thus, if the PO has e p > e c /2, orbits with e 0 < 2e p -e c (e.g. a circular orbit) will crash on Mars surface in less than T s /2 ; this is essentially what Knežević and Milani (1998) found, while studying the crash of the Apollo 16 subsatellite on the Moon. This is exemplified in figure 9, for two quasi-circular low-altitude orbits around Mars at inclinations i 0 = 90 • and i 0 = 60 • , both of which crash in less than T s (i 0 )/2. Note that the slow period T s (i 0 )/2 depends strongly on i 0 , becoming very large at i c ≈ 63.4

• . Yet, this period will be only 25 Martian days at i 0 = 90 • and at most 100 day ⊗ as soon as |i 0 -i c | > 3 • . On the other hand, the PO itself will keep its eccentricity nearly constant and will never crash. Another example is shown in figure 10 where a slowly evolving orbit close to the critical inclination i c starting at low eccentricity nevertheless crashes after 79 day ⊗ (about 3 times the orbit arc used in Prony's FA), as the PO is very close to the collision limit. Choosing a circular orbit or the FO of the 'J 2 &J 3 problem' as initial condition, as is often done for preliminary design, the orbit would again crash, showing that these approximations should always be complemented by other methods including higher-degree terms of the gravitational potential at low altitudes, as also advocated in [START_REF] Lara | Precise Analytical Computation of Frozen-Eccentricity, Low Earth Orbits in a Tesseral Potential[END_REF]. The true PO found by our algorithm is indeed stable and safe from collision.

For orbits starting with opposite sign for g 0 with re- ℑ(α 0 ) ± α p FO ℑ(α 0 ) PO e 0 sin(g 0 )

Figure 11: Dynamical fate map, showing the survivability of lowaltitude Martian orbits, in the i 0 -e 0 sin g 0 plane. Blue circles are the POs; the dashed line represents the orbits mean position e F sin g F (FOs). These orbits are long-lived, as eccentricity excursions are within the turquoise band (±αp). Orbits in the white zones are also long-lived, with eccentricity excursions of order 2 δe F around the PO. The magenta zones contain initial conditions that collide with Mars in less than Ts/2 + Trot. In the dotted zones, orbits will collide in less than Ts/2 while, in the red zones, orbits would collide in less than Trot. In the grey zones, orbits are already above ec. spect to g p , the situation is worse, as |α s | = |e 0 + e p |, i.e. e 0 is pushed up to 2e p + e 0 and the satellite crashes if this value is larger than e c . Hence, orbits with e 0 > e c -2e p will crash on Mars and, for e p ≥ e c /2, all orbits with that (wrong) sign of g 0 -even with e 0 ≈ 0 -will also crash. Moreover, if e p reaches e c , any orbit will crash in less than T s /2 whatever its initial eccentricity e 0 or the sign of its initial pericenter argument g 0 .

From the above analysis it is evident that the PO position conditions the existence of non-collisional orbits at low altitudes, as well as the extent to which active control needs to be used, in order to control eccentricity and altitude variations.

Dynamical fate maps

Using the results presented in the previous sections, we can build 'dynamical fate maps' that immediately display the survivability of orbits around POs. An example is given in figure 11, as a projection on the (i 0 , e 0 sin g 0 ) plane, for a = R ⊗ + 60 km.

Let us describe all orbital regions, shown on this map by different color/texture. We know that a PO at given i 0 is stable and does not lead to collision as, in general, the periodic terms have amplitudes of order α p 10 -3 (the turquoise band in figure 11). To fix ideas, consider g p > 0, as is the case in figure 11 for 0 < i 0 < 13 • . If we move to 'higher' positions e c > e 0 > e p (white zone above the PO), the orbit will oscillate between e 0 and 2e p -e 0 (with the latter negative for large e 0 ) but will never reach the (lower) collision eccentricity -e c , so the orbit will live forever with eccentricity variations |e 0 -e p |. If we now move to 'lower' positions e 0 sin g 0 , including those for which g 0 is negative (white zone below the PO), the orbits will still live forever, but will oscillate in e with an amplitude |α s | ≈ |e 0 sin(g 0 )e p |. This amplitude increases when moving lower, until we reach position -e c +2e p +α p , where the combined effect of proper and forced terms will cause the orbit to crash after at most T s /2 + T rot (magenta zones in figure 11). This zone has a height in eccentricity α p and, as soon as we reach -e c +2e p , orbits will collide in less than T s /2, simply due to the slow proper motion (dotted zones in figure 11). Below e 0 = -e c + α p or above e c -α p , the small periodic terms are enough to cause the orbit to crash in at most one rotation of Mars T rot .

The map provides a rather complete picture of the dynamical fate, for all initial conditions at a given altitude. It explains both collision zones observed for e 0 < e c in figure 5, i.e. for 0 < i 0 < 13 • and 75 • < i 0 < 105 • . It also shows that, for 61 • < i 0 < 65 • , non-collisional orbits do not exist, as these POs would have e p > e c .

It is interesting to see how the map changes with altitude. There is no technical difficulty in building this map for any altitude. Maps of the entire space of initial conditions a 0 , i 0 , e 0 can be obtained efficiently using our method, but it is difficult to visualize the results in 3-d. We thus display a 2-d projection of these results in the (a 0 , e 0 sin g 0 ) space, for selected inclinations i 0 in figure 12. Note that the initial conditions for the POs are expected to be smooth functions of a 0 , as for i 0 . In fact, this implies that the whole 3-d map can be computed by continuation, starting only from a single guess, locating a PO at some altitude and i 0 , then using this as initial guess for an adjacent value of i 0 , etc. When the inclination range is exhausted, we can start incrementing the altitude, repeating the above procedure. The displacement of FOs (i.e. mean positions of POs) and of the initial conditions of POs with altitude is nontrivial. While the size of the fast terms α p is decreasing monotonically with altitude, the FOs displacement is more complicated; this is not so surprising as in the axisymmetric case the FO location can be found by numerically solving a strongly non-linear equation, whose coefficients depend on altitude (see e.g. Lara et al. (2009)). Still, this shift is small at low altitudes, while e c is proportional to a 0 -R ⊗ . Thus, there is always an altitude for which e p ≥ e c and no orbit can survive below. A dynamical fate map for polar orbits is displayed in figure 12 as a function of altitude. As one can see, polar Martian orbits are safe down to a 0 = R ⊗ + 33 km, with e ≈ 0.008. This orbit resembles the one of the Global Mars Surveyor (GMS), but at a much lower altitude. Note that, for polar orbits, e p varies very little with altitude, so it remains roughly the same also at the GMS altitude (372 km).

Near-frozen, quasi-circular orbits at very low altitudes ( 10 km) are possible, but only for selected inclinations. As seen in figure 12, this is the case for i 0 ≈ 13 • and i 0 ≈ 66 • . The dynamical fate map for i 0 = 66 • is shown in figure 13. In principle, non-crashing orbits down to 6 km can be found, but the Martian topography would certainly not allow it.

Conclusions -Discussion

In this paper we presented a detailed study of low Martian orbiters (a -R ⊗ < 500 km), applying Prony's Frequency Analysis and an iterative filtering algorithm that allowed effective computation of families of periodic orbits that dominate the secular dynamics of nearly circular satellite orbits. Our main goal was to present an effective way of describing eccentricity variations and determining the dynamical fate of orbits, such that it could be useful for preliminary orbit design. As high-eccentricity orbits cannot have near-constant altitude, our primary focus was on the 2-d.o.f. equivalent of frozen orbits, namely nearlyfrozen periodic orbits.

By performing a series of tests, we found the 'minimal' model for the Martian gravity field to be of degree l max = 20 and studied both the full model (20DM) and its axisymmetric counterpart (20ZM). In contrast to the Moon, we find that the Martian gravity field can be well approximated by an axisymmetric model, even for satellites at altitudes 60 km, with the exception maybe of nearequatorial orbits (i ≤ 13 • ). Nevertheless, the true secular motion is described by two vectors (the complex u(t) and v(t) in our notation), although weakly coupled in most cases, and an accurate representation of both should be given.

In the absence of exact analytical solutions for the highorder and degree 2-d.o.f. problem, Prony's frequency analysis is a very effective method of accurately computing a quasi-periodic fit of the secular orbit. The separation of the two main frequencies guarantees that low-order linear combinations (resonances) will not 'contaminate' the signal and the motion will be regular. Prony's method has some advantages over other spectral methods, namely it requires only a very short orbit arc to produce an accurate determination of the fundamental modes (frequencies and amplitudes). Moreover, under the same assumption of frequency separation, it can be easily complemented by an iterative filtering method, which allows efficient determination (2-3 iterations) of periodic orbits, in the vicinity of any initial condition. Note that, apart from locating the POs, it also provides a very precise decomposition of these orbits, which contains only low-amplitude, 'fast', periodic variations around the mean position. The constructive theory of [START_REF] Wiesel | A Theory of Low Eccentricity Earth Satellite Motion[END_REF] explains the frequency decomposition of POs and provides ground for the convergence of our iterative algorithm. Unlike differential corrections, our method unfortunately does not provide direct indication of linear stability for the POs. On the other hand, the correction steps are independent of the complexity of the gravity model considered.

A collection of POs for different altitudes and inclinations, projected on the (e cos g, e sin g) plane, is given in figure 14. Because PO positions are smooth functions of the initial conditions, one can compute a full 3-d map of POs (in a, e and i) by continuation, starting from a first guess and smoothly varying the parameters of the family of POs. While the overall computational burden is certainly higher than e.g. solving a nonlinear equation of high degree to find the FOs of an axisymmetric model, the gain in precision, but mostly in acquiring also the spectral decomposition of the true secular orbit, is considerable. Using the information on relative amplitudes, one can construct and fully describe 'dynamical fate maps' for satellite orbits at any given altitude and inclination. These maps can be seen as a generalization of the classical inclinationeccentricity diagrams and are a potentially useful tool for preliminary mission design. Figure 12, for example, confirms that the Global Mars Surveyor orbit (i = 93 • , e = 0.008, minimum altitude 372 km) keeps very close to the corresponding PO and has eccentricity variations of only 10 -3 -smaller than at most other inclinations for the same altitude -and altitude variations of ≈ 60 km. The Exomars Trace Gas Orbiter (i = 74 • , e = 0, altitude 400 km) and Mars Odyssey (i = 93 • , e = 0, altitude 400 km) require active control to suppress altitude variations, as they are located at δe F ≈ 0.008 from the PO, and hence have eccentricity variations 2 δe F ≈ 0.016 and altitude variations of ≈ 120 km; still, the collision eccentricity is ≈ 0.1 at 400 km altitude and so these satellites are safe. However, figure 11 suggests that the same strategy (i.e. choosing e = 0) would become inefficient at lower altitudes and, in fact, for i = 74 • or 93 • and initial altitudes ≤ 70 km, a circular orbit would swiftly crash on Mars.

In conclusion, we believe that our method can be an effective tool for preliminary trajectory design. We plan to explore it further -considering other solar system objects as targets -and report in a forthcoming paper.

Acknowledgements

We are grateful to the anonymous reviewers for their thorough reviews that helped us improve our paper. Many thanks go to S. Tzirti for assisting us with numerical integrations used in this study, even long after having left the field. . With decreasing altitude, the PO's mean position (FO, crosses) can move to higher (for i 0 = 7 • , 61 • or 90 • ) or lower eccentricities (i 0 = 66 • ). At lower altitudes, the periodic terms increase in amplitude and the projection of the orbit becomes more and more involved; FA can accurately describe this behavior. For i 0 = 61 • and 50 km altitude the satellite would crash on Mars.

  al. (2009); Tzirti et al. (2010); Tzirti et al. (2014); Noullez et al. (2015)

Figure 1 :

 1 Figure 1: Normalized gravitational coefficients J l,m ≡ (C l,m 2 + S l,m 2 ) 1/2 for the Moon (left) and Mars (right) shown at the same scale, as a function of degree l. Different orders m are shown with different symbols and colors (red, violet, magenta). Zonal coefficients m = 0 are shown as blue circles, empty • for even degree l and filled • for odd degree. The thick dashed line is the degree mean amplitude h P l m=0 J l,m 2 /(2l + 1) i 1/2as a function of degree l, and decays approximately as l -2 (Kaula's rule, see[START_REF] Kaula | Theory of Satellite Geodesy[END_REF]).

a

  = R ⊗ + 60 km, g 0 = -90°, h 0 = 0, e 0 = 0.00252, 20ZM/20DM 20DM 20ZM ν 1 (J 2 )

a

  = R ⊗ + 60 km, g 0 = -90°, h 0 = 0, e 0 = 0.00252, 20ZM/20DM 20DM 20ZM ν s (J 2 )

aa

  Figure6: One-dimensional cuts in the map of initial conditions for the FA amplitudes, on the left along i 0 (horizontally) for a constant initial e 0 = 0.00873 (shown as dashed thin horizontal line on the left), and on the right along e 0 (vertically) for i 0 = 90 • , both for fixed initial orientation angles g 0 = -90 • and h 0 = 0. The three curves in the two graphs are the slow mode amplitude |αs| (blue circles), the constant amplitude e F = |α 0 | (black squares, empty if g F = -90 • , filled if g F = +90 • ) and the total periodic amplitude αp (red diamonds). Around i 0 ∼ 63 • , orbits collide with Mars in less than 30 day ⊗ and cannot be analyzed. The dashed thick violet lines on the left and right correspond to the collision eccentricity ec = 0.0174 at that altitude. Orbits will collide with Mars surface (crossed or hatched zones on the right) if |α 0 | ± αp (grey zone on the right), added to ±|αs| (thin blue lines on the right) reaches the collision eccentricity ec, which can occur at either high or low initial eccentricity e 0 .

3 Figure 10 :

 310 Figure10: Orbital evolution (thin black lines and circles) for i 0 = 61 • and e 0 = 0.002; the orbit collides with Mars after 79 day ⊗ (violet dashed circle denotes the collision line). Another orbit, started at e 0 = 0.007 (the FO of the 'J 2 &J 3 problem') also crashes after 88 day ⊗ (dotted magenta curve). The PO obtained after 3 iterations of our filtering algorithm (red line) will never crash on the Martian surface, and stays around the 'nearly-frozen' parameters g = -89 • , e = 0.0141, having excursions in eccentricity ≈ ±0.0005.

a

  = R ⊗ + 60 km, h 0 = 0, 20DM(3)

  Displacement of FOs (i.e. mean positions of POs) locations as a function of the altitude a 0 -R ⊗ for different initial inclinations i 0 = 7 • , i 0 = 61 • , i 0 = 66 • and i 0 = 90 • . Orbits in the shaded grey zones exceed ec (thick, dashed violet line). The POs (blue circles) and the dynamical fate zones are plotted for i 0 = 90 • (the bands of width ±αp are invisible at this scale). The PO eccentricity reaches ep = ec at about 33 km altitude for polar orbits.

  0 = 66 °FO i 0 = 66°P O i 0 = 66°Figure 13: Dynamical fate map for i 0 = 66 • . The PO is located at nearly zero eccentricity. Symbols and zones color codes are the same as in figure11. Orbits could survive down to 6 km altitude.

  0 = 0, i 0 = 7°, 20DM(1)a = R ⊗ + 50 km a = R ⊗ + 200 km a = R ⊗ + 300 km a = R ⊗ + 500 km 0 = 0, i 0 = 66°, 20DM(1) a = R ⊗ + 50 km a = R ⊗ + 100 km a = R ⊗ + 200 km a = R ⊗ + 300 km a = R ⊗ + 400 km a = R ⊗ + 500 km -0 = 0, i 0 = 61°, 20DM(1) a = R ⊗ + 50 km a = R ⊗ + 200 km a = R ⊗ + 300 km a = R ⊗ + 400 km a = R ⊗ + 500 km 0 = 0, i 0 = 90°, 20DM(1) a = R ⊗ + 50 km a = R ⊗ + 300 km a = R ⊗ + 500 km

Figure 14 :

 14 Figure14: Periodic Orbits for different values of a 0 -R ⊗ < 500 km and i 0 : 7 • (top left), 61 • (bottom left), 66 • (top right) and 90 • (bottom right). With decreasing altitude, the PO's mean position (FO, crosses) can move to higher (for i 0 = 7 • , 61 • or 90 • ) or lower eccentricities (i 0 = 66 • ). At lower altitudes, the periodic terms increase in amplitude and the projection of the orbit becomes more and more involved; FA can accurately describe this behavior. For i 0 = 61 • and 50 km altitude the satellite would crash on Mars.

  Frozen or Sun-synchronous or Repeating Ground Track) orbits around Mars are rather rare (but see the recent papers of Liu et al. (2010), Liu et al. (2011), Liu et al. (2012) or Zhou et al. (2012)), as it is only recently that accurate models of Mars gravity field have been obtained (see Konopliv et al. (2016), Genova et al. (

Table 1 :

 1 Main physical parameters for the Earth, Moon and Mars.

		Earth	Moon	Mars
	Mass M [kg]	5.972 × 10 24	73.48 × 10 21	641.9 × 10 21
	Radius R [km]	6378	1738	3397
	Oblateness J2	1.083 × 10 -3	204 × 10 -6	1.956 × 10 -3
	J3 -2.532 × 10 -6 9.999 × 10 -6 31.94 × 10 -6
	Trot [day]	1	27.32	1.026
	Satellite at body surface a = R
	eF(J2 &J3; R, 90 • )	0.0012	-0.0245	-0.0082
	Trev(R) [day]	0.0587	0.0752	0.0696
	Ts(R) [day]	18.06	123	11.86
	Trot/Trev(R)	17	363.2	14.75
	Ts(R)/Trot	18.06	4.5	11.56
	Satellite at 100 km altitude a = R + 100 km
	ec(a)	0.0154	0.0544	0.0286
	Trev(a) [day]	0.06	0.0818	0.0727
	Ts(a) [day]	19.08	149.5	13.12
	Trot/Trev(a)	16.6	334	14.12
	Ts(a)/Trot	19.08	5.7	12.79

From now on, we will express times in units of Martian sideral day ⊗ ≡ Trot = 1.025957 Earth days, and frequencies in units of day ⊗ -1 .

In this section, we restrict ourselves to cases with g 0 = ±90 •


Gang Zheng
email: gang.zheng@imag.fr

Arnaud Tonnelier
email: arnaud.tonnelier@inrialpes.fr

Chaotic solutions in the quadratic integrate-and-fire neuron with adaptation

Keywords: QIF, Chaos, Poincaré map, Lyapunov exponent

The quadratic integrate-and-fire (QIF) model with adaptation is commonly used as an elementary neuronal model that reproduces the main characteristics of real neurons. In this paper, we introduce a QIF neuron with a nonlinear adaptive current. This model reproduces the neuron-computational features of real neurons and is analytically tractable. It is shown that under a constant current input chaotic firing is possible. In contrast to previous study the neuron is not sinusoidally forced. We show that the spike-triggered adaptation is a key parameter to understand how chaos is generated.

Introduction

Neurons transform incoming stimuli into train of all-or-one electric events known as spikes. Precise spike times are believed to play a fundamental role in the encoding of information by neural systems. In this view, the spike time reliability is a key to understand the basis of the neural code [START_REF] Mainen | Reliability of spike timing in neocortical neurons[END_REF]. It has been shown that the responses of neurons to time varying stimuli have high reliability with a precision of 1 ms or less [START_REF] Mainen | Reliability of spike timing in neocortical neurons[END_REF], [de Ruyter van Steveninck et al., 1997], [START_REF] Berry | The structure and precision of retinal spike trains[END_REF]. On the other hand, constant stimuli elicits spike trains with high variability. It has been primarily suggested that such variability in discharge times is due to accumulating noise. However, it could be due as well to a chaotic behavior [START_REF] Hayashi | Chaotic behavior in the onchidium giant neuron under sinusoidal forcing[END_REF], [START_REF] Chacron | Chaotic firing in the sinusoidally forced leaky integrate-and-fire model with threshold fatigue[END_REF]. In fact, the responses of a neuron to given stimuli over repeated trials could be different because of intrinsic noise but also because of the chaotic behavior of the neuronal dynamics.

The reproducibility of spikes has been largely investigated in theoretical studies. Simplified phenomenological neuron models of integrate-and-fire type are commonly used to explore the dynamical properties of neurons and particularly the robustness of neuron response [Tiesinga, 2002], [START_REF] Brette | Reliability of spike timing is a general property of spiking model neurons[END_REF]]. Chaotic solutions have been investigated for the leaky integrate-and-fire neuron [Coombes, 1999], [START_REF] Chacron | Chaotic firing in the sinusoidally forced leaky integrate-and-fire model with threshold fatigue[END_REF]. A time-varying injected currents combined with a periodic modulation of parameters are the minimal requirements to obtain chaotic solutions.

The leaky integrate-and-fire neuron has some obvious limitations and an extension to quadratic integrate-and-fire (QIF) neurons allows to describe the nonlinear spike-generating currents of real neurons. In this paper we show that adaptive QIF neuron can exhibit chaotic solutions under constant input current. This result corroborates the non-reliability of spike times observed experimentally when a constant current is injected into a neuron [START_REF] Mainen | Reliability of spike timing in neocortical neurons[END_REF], [de Ruyter van Steveninck et al., 1997].

The paper is organized as follows. In section 2 the model is presented. We derive in section 3 the expression of the Poincaré map associated with the firing of a spike and use the Marotto theorem in section 4 to show that the adaptive QIF can exhibit chaotic dynamics. Numerical simulations are made in section 5 to support our analysis and to highlight the dynamical behavior of the model.

The model

We consider an adaptive quadratic integrate-and-fire model given by the following differential equations:

dx dt = x 2 + a -y (1)
dy dt = b -y τ (x) (2)
and the spike rules:

x(t + ) = q, if x(t) = h (3) y(t + ) = cy(t) + p, if x(t) = h (4)
where x is the membrane potential and y a phenomenological variable that represents adaptation. Parameters a, h and q are respectively the constant input current, the peak of a spike (see below) and the reset value. Parameters b, c and p describe the adaptive current and we assume c ≥ 0. Note that y can be seen as a recovery variable. The function τ (x) models the voltage-dependence of the adaptive variable time-constant. For a fixed voltage value x, the variable y approaches the equilibrium b with a time constant τ (x). In the following we consider τ (x) = τ /x where τ is a positive constant. This choice approximates the voltage dependence of the time constant of the adaptive variable commonly used in detailed models: for large value of x the integration time of y is shorter. Therefore the time evolution of the adaptive current is given by

dy dt = x(b -2y)/τ
where the factor 2 is used for convenience (see below) and can be removed (change of variables).

Because of the quadratic term in (1) the membrane potential can escape to infinity in finite time. We thus introduce h the cutoff of a spike that defines the firing of a spike. After a spike has been triggered the membrane potential restarts at a reset value q and the variable y is multiplied by a constant c and is increased by an amount p which accounts for spike-triggered adaptation. The differential equations (1)-(2) model the subthreshold dynamics and the spike initiation. When x attains the peak, the impulsive effect occurs according to (3)-(4). The time evolution is illustrated in Figure 1 for parameter values where the model exhibits regular spiking. Simulations (not shown) indicate that classical neuron-computational features of neurons [Izhikevich, 2004] including bistability, resonator and bursting (see later) can be reproduced by this model. For simplicity but without loss of generality, we consider τ = 1 in the following analysis.

Poincaré map

Due to the non smooth nature of integrate-and-fire models, the standard theory of dynamical systems does not apply directly. This problem was overcome in previous studies using leaky integrate-and-fire model where a closed form for the map of successive firing times can be analytically derived. However, because of the nonlinearity of the subthreshold dynamics of the QIF it is no longer possible here to derive an analytical expression for this map. In order to analyze the behavior of system (1)-( 4), we introduce the following Poincaré section: S = {(x, y) | x = h} . It is easy to show that the differential system (1)-( 2) is conservative. In fact, if we define the energy of the system as follows:

E(x, y) = 1 2 y 2 -y(a + x 2 ) + b 2 x 2 (5)
we have dE/dt = 0 along the trajectory of the system provided that no firing occurs. Let us now consider a trajectory that reaches the section S (the neuron fires) at point P k = (h, y k ), then according to the spike rules ( 3)-( 4), the trajectory jumps to the point Pk = (q, cy k + p) due to the impulsive effect (see Figure 2). After that the trajectory evolves according to (1)-( 2), and it could reach again the section S at a point P k+1 = (h, y k+1 ). Since Pk and P k+1 are both located on the trajectory of system (1)-( 2) where no impulsive effect occurs, these two points satisfy relation ( 5). Consequently, we obtain:

1 2 y 2 k+1 -a + h 2 y k+1 + b 2 h 2 = 1 2 (cy k + p) 2 -(q 2 + a) (cy k + p) + b 2 q 2
and some straightforward calculations yield:

y k+1 -a -h 2 2 = cy k + p -a -q 2 2 + 2a + h 2 + q 2 -b h 2 -q 2
that defines an implicit map for the adaptive variable at the firing times. Assuming that ∀k, y k < a + h 2 and setting L = (2a

+ h 2 + q 2 -b)(h 2 -q 2 ), H = a + h 2 , Q = p -a -q 2 ,
then we can define the map f

y k+1 = f (y k ) = H -(cy k + Q) 2 + L (6)
that gives the successive values of the adaptive variable at the firing times.

Chaotic behavior

In this section, with the help of Marotto's Theorem [Marotto, 1978], [Marotto, 2004] we prove that the adaptive QIF neuron ( 1)-( 4) can exhibit chaotic behavior for specific values of parameters. The main idea of Marotto's theorem is to seek a snap-back repeller. Marotto showed that the presence of a snap-back repeller is a sufficient criterion for the existence of chaos. Let us first recall the basic Marotto's Theorem for one-dimensional system:

Theorem 1 [Marotto, 1978] Let y * be a fixed point of the map f . This fixed point is said to be a snap-back repeller if the following conditions are satisfied: A snap-back repeller arises when the fixed point is repelling while there is a trajectory starting in the repelling neighborhood of the fixed point that goes back to the fixed point, i.e. the repelling fixed point has an associated homoclinic orbit. Note that the existence of a snap-back repeller for one-dimensional system is closely related to the existence of a periodic point with period 3. In fact, it can be shown that the existence of a snap-back repeller of f is equivalent to the existence of a point of period 3 for the map f n for some positive integer n [START_REF] Li | Period three implies chaos[END_REF].

(i) f is differentiable in a neighborhood B r (y * ) of
Using Theorem 1 we will prove the following Proposition 1 For appropriate values of parameters the map f is chaotic in the sense of Marotto.

In the remainder of this section we prove that the map f fulfills the conditions of Marotto theorem. In order to identify a snap-back repeller for the one dimensional map (6), we derive the eigenvalue of the system and show that, for values of the state variable close to the fixed point, the corresponding eigenvalues lie outside the unit circle, whereas for values of the state variable that are sufficiently far from the fixed point, the eigenvalues are within the unit circle.

Let us firstly check whether condition (i) in Theorem 1 is satisfied. From (6) it is easy to show that the map has a fixed point y * = f (y * ) if the following inequality is fulfilled:

(cQ + H) 2 -Q 2 -H 2 + L c 2 -1 ≥ 0 (7)
which yields the following two fixed points:

y * = ± (cQ + H) 2 -(Q 2 -H 2 + L) (c 2 -1) -(cQ + H) c 2 -1 (8)
Moreover, according to (6) we have

Df (y) = -c cy + Q (cy + Q) 2 + L
where Df represents the derivative of f . The function Df is decreasing if

L ≥ 0 (9)
since c > 0. In addition, by defining y z = -Q c it is easy to check that f (y) is increasing if y < y z and is decreasing when y > y z .

Let y c be the solution of the equality: Df

By construction, for all y ∈ B r (y * ), Df (y) is always strictly larger than one in absolute values provided that y * > y c > y z . Therefore, condition (i) of Theorem 1 is satisfied.

On the other hand, in order to fulfill condition (ii) of Theorem 1, we seek a point y m = y * in the neighborhood B r (y * ), such that y m-1 = f (y m ), . . ., y 1 = f (y 2 ) = f m-1 (y m ), and y * = f (y 1 ) = f m (y m ), with some positive integer m, and Df (y k ) = 0, for 1 ≤ k ≤ m.

Let us consider the following fixed point

y * = (cQ + H) 2 -(Q 2 -H 2 + L) (c 2 -1) -(cQ + H) c 2 -1 (12)
and its neighborhood B r (y * ).

According to map (6), we have y * = f (y 1 ) = H -(cy 1 + Q) 2 + L which yields two possible solutions:

y 1,1 = (y * -H) 2 -L -Q /c and y 1,2 = -(y * -H) 2 -L -Q /c
. The first solution is identically equal to its corresponding fixed point, hence the possibility of existence of a snap-back repeller y m ∈ B r (y * ) is relative to the second one, noted as y 1 . Obviously we have y 1 < y * .

Assume that the following conditions are satisfied:

H > f (y z ) = y A > y * f (y A ) = y B < y 1 (13)
Since f (y A ) = y B and f (y * ) = y * , and f (y) is a decreasing continuous function over [y * , y A ], therefore there always exists a point y 2 ∈]y * , y A ], such that y 1 = f (y 2 ). If such a y 2 is located in B r (y * ), we have m = 2. Otherwise, consider the interval ]y 1 , y z [, since y * = f (y 1 ) and y A = f (y z ) , and f (y) is increasing over [y 1 , y z ], it is always possible to find a y 3 ∈ [y 1 , y z ] such that y 2 = f (y 3 ).

It should be noted that there exists another possible y 3 satisfying y 2 = f (y 3 ), which is located over [y z , y * ]. Notice that y * = f (y * ), y A = f (y z ), and f (y) is decreasing over [y z , y * ], thus we can find the second y 3 ∈ [y z , y * ] such that y 2 = f (y 3 ). Again if y 3 ∈ B r (y * ), then we found m = 3. Otherwise, let us firstly consider In order to prove the convergence of the above sequence to the fixed point y * , let us analyze the right side of the sequence (i.e., even subscript), which in fact can be written as follows:

y 3 ∈ [y 1 , y z ]. With f (y 2 ) = y 1 , f (y * ) = y * , we get y 4 ∈ [y * , y 2 ] such that f (y 4 ) = y 3 . If this y 4 is not yet located in B r (y * ), according to f (y 3 ) = y 2 , f (y z ) = y A , f (y * ) = y * , since f (y) is
y 2k = f (f (y 2k+2 )) = H - c H -(cy 2k+2 + Q) 2 + L + Q 2 + L with k ≥ 1. When k → ∞,
we have y 2k = y 2k+2 , noted as y lim . Obviously y lim one solution of equation: y lim = f 2 (y lim ), which in fact has two solutions under the constrain y * < y 2k+2 < y 2k < H. It is easy to check that equilibrium points of ( 8) are these two solutions. Since our analysis is focused only on the equilibrium point of ( 12), thus we have y lim = y * . It should be noticed that, for the left side of the sequence (odd subscript), we can obtain the same result through similar process, which implies the sequence converges towards the fixed y * via both two sides. Consequently, we can always find a point y m ∈ B r (y * ), such that f m (y m ) = y * with some positive integer m.

With the same procedure, it is easy to adapt the above proof for the second possible value of y 3 ∈ [y z , y * ].

Finally, since Df (y k ) = -c

cy k +Q √ (cy k +Q) 2 +L
, hence if

y k = - Q c ( 14
)
we have Df m (y k ) = 0. It is worth noting that if there exist some points such that y i = -Q c , then it corresponds to a periodic orbit of ( 6). In summary, provided that the technical conditions: ( 7), ( 9)-( 10), ( 13)-( 14) are fulfilled, it is possible to determine a neighborhood B r (y * ) around the fixed point ( 12) with r defined in (11) such that conditions (i) and (ii) are satisfied. Therefore, Theorem 1 holds and this ends the proof of Proposition 1.

Numerical simulations

Besides the theoretical analysis, numerical simulations are made in order to show that system (1)-( 4) can exhibit chaotic dynamics. In particular, we show that there exists parameter values for which the technical assumptions used in our proof are satisfied.

We consider the following QIF neuron with adaptation:

dx dt = x 2 + 6 -y (15) dy dt = x(2 -2y)
and the spike rules

x(t + ) = 10, if x(t) = 20 (16) y(t + ) = cy(t) -0.2, if x(t) = 20
i.e. we choose the following set of parameter values: a = 6, b = 2, τ = 1, p = -0.2, q = 10 and h = 20. We explore the influence of the parameter c on the asymptotic behavior of the model and we investigate the existence of adaptation-induced chaotic dynamics using different values for c.

For c = 13.8, we have L = 153000, H = 406 and Q = -106.2 that give the corresponding Poincaré map:

y k+1 = 406 -(13.8y k -106.2) 2 + 153000 (17)
Numerically we check that inequality ( 7) is satisfied and straightforward calculations yield the following values for the variables introduced in the proof:

y * ≈ 11.4434, y c ≈ 9.7957 y z ≈ 7.6957, y 1 ≈ 3.9479, y A ≈ 14.8479, y B ≈ 2.5873 It is easy to check that conditions ( 9)-( 10) and ( 13)-( 14) are all satisfied. Moreover we choose r = 1.5 < y * -y c that defines the neighborhood B r (y * ) = ]9.9434, 12.9434[ (see Figure 3A). Since we have Df (9.9434) = -1.0909 and Df (12.9434) = -2.5123, thus for all y k ∈ B r (y * ), condition (i) of Theorem 1 is fulfilled.

It is easy to check that, for a y 4 ≈ 12.6150 ∈ B r (y * ) and y 4 = y * , we have y 3 = f (y 4 ) ≈ 9.0005, y 2 = f (y 3 ) ≈ 14.4336, y 1 = f (y 2 ) ≈ 3.9479 and y * = f (y 1 ) (see Figure 3B). Moreover Df 4 (y 4 ) ≈ -8.6461 = 0, hence condition (ii) of Theorem 1 is satisfied. Consequently, the point y * is a snap-back repeller and the map (6) exhibits chaotic behavior.

Besides the existence of chaos, the scenario leading to such dynamics as c varied is illustrated in Figure 4. The phase portrait of system (15)-( 16) is shown together with the membrane potential as a function of time. When c = 10, the system tends towards an asymptotically stable periodic orbit (See Figure 4A,B). In Figure 4C,D, for c = 13.9 a period 3 solution occurs. A chaotic behavior can be generated by choosing c = 13.8 as shown in Figure 4E,F.

To gain deeper understanding of the chaotic dynamics, the bifurcation diagram obtained while varying parameter c is shown in Figure 5A, while Figure 5B shows the leading Lyapunov exponent of the system. The lines in the bifurcation diagram correspond to periodic spiking and the thick parts are associated with regimes where the adaptive variable takes either a countable but very large number of values or an uncountable value. It is not clear that in the bifurcation diagram the regimes where the adaptive variable takes a very large number of values are associated with periodic orbits with large period or chaotic regimes.

One of the usual ways to test for chaos is to calculate system's largest lyapunov exponent. The Lyapunov exponent of a dynamical system is a quantity that characterizes the rate of separation of infinitesimally close trajectories. The general idea is to follow two nearby orbits and to calculate their average logarithmic rate of separation, and a positive Lyapunov exponent indicates chaos.

According to expression (17), it is easy to calculate the largest lyapunov exponent and the result is shown Figure 5B. The range of positive exponent suggests chaotic regimes within these regions. Moreover the bifurcation diagram reveals that the system exhibits transition to chaos through period doubling bifurcation cascade.

Conclusion

Periodically forced integrate-and-fire neurons have been proposed to produce chaotic firing. It is expected that chaos in neuronal models can be generated through an applied constant current in agreement with the experimentally observed neural spike trains variability in this condition. This paper presents a simple nonlinear integrate-and-fire neuron which can exhibit a chaotic behavior under a constant current. The model is the quadratic integrate-and-fire neuron with a nonlinear adaptive current that takes into account the voltage dependence of the adaptive-variable time constant. The situation significantly differs from those in previous studies that report chaos in integrate-and-fire models for several reasons: (i) there is no periodic forcing (ii) threshold or resetting voltage is not periodically modulated and (iii) the integrate-and-fire model is nonlinear. We proved mathematically the existence of chaotic solutions provided that parameters of the adaptive current are suitably chosen. Numerical simulations are given including the largest Lyapunov exponent and the bifurcation diagram.

A further development of this work is the study of entrainment in sinusoidally driven adaptive QIF neuron. It is expected that the model exhibits specific dynamics in response to periodic forcing that do not exist in LIF models. Further research is needed to clarified this point.

Figure 1 :

 1 Figure 1: Illustrative plots of the time evolution for A: the membrane voltage x(t), and for B: the adaptive variable y(t) in the regular spiking regime. Parameters are a = 2, τ = 15, b = 1, q = 1, h = 8, c = 1 and p = 0.25.

Figure 2 :

 2 Figure2: A piece of trajectory containing a jump (spike) in the phase plane (x, y). The vertical dashed lines are the reset line x = q and the Poincaré section that corresponds to the threshold line x = h.

  the fixed point y * with radius r > 0 and the eigenvalue of Df (y) is strictly larger than one in absolute value for all y ∈ B r (y * ). (ii) There exists a point y m ∈ B r (y * ), with y m = y * , such that for some positive integer m, f m (y m ) = y * , and f m is differentiable at y m with Df m (y m ) = 0.

  (y c ) = -1. Since lim y→-∞ Df = c and lim y→+∞ Df = -c, a necessary and sufficient condition for the existence of y c is c > 1. (10) Moreover since Df (y c ) < Df (y z ) and Df (y) is decreasing when y > y z , we have y c > y z . Let us now consider a neighborhood of y * , noted as B r (y * ), with 0 < r ≤ min{|y * -y c | , H -y * }

  increasing over [y 3 , y z ] but decreasing over [y z , y * ], hence it always exists a y 5 ∈ [y z , y * ] such that f (y 5 ) = y 4 . By induction, if the found points are not yet entered into B r (y * ), the above procedure generates the sequence: y B < y 1 < y 3 < y z < y 5 < • • • < y * < • • • < y 4 < y 2 < y A < H, which tends towards the equilibrium point y * .

Figure 3 :

 3 Figure 3: Snap-back repeller of the QIF neuron with adaptation. The numerical investigations for the existence of a snap-back repeller are done in A and B, respectively. A: The derivative of the Poincaré map (17). The successive iterations y k are represented in the neighborhood of the fixed point y * . Numerically, we check that the absolute value of the derivative is always greater than 1 for all y k . B: Location of the iterations with respect to the neighborhood of the fixed point y * . We find y 4 ∈ B r (y * ) and y 4 = y * together with f 4 (y 4 ) = y * .

Figure 4 :

 4 Figure 4: Transition to chaos through a cycle of period 3 in the QIF neuron with adaptation. The phase portraits are given in the left column and the x dynamics towards time are depicted in the right column. A, B: Asymptotically orbit for system with c = 10. Simulation was made with initial conditions: x 0 = 5, y 0 = 15. C, D: Period 3 orbit obtained for c = 13.9. Initial conditions are x 0 = 15, y 0 = 15. E, F: Chaotic orbit for system with c = 13.8 with initial conditions: x 0 = 10, y 0 = 10.