N Delsate

• P Robutel

• A Lemaître

• T Carletti

Frozen Orbits at high eccentricity and inclination: Application to Mercury orbiter

Keywords: Methods: analytical study, Stability, Long-term evolution, Kozai resonances, Frozen Orbit equilibria

We hereby study the stability of a massless probe orbiting around an oblate central body (planet or planetary satellite) perturbed by a third body, assumed to lay in the equatorial plane (Sun or Jupiter for example) using a Hamiltonian formalism.

We are able to determine, in the parameters space, the location of the frozen orbits, namely orbits whose orbital elements remain constant on average, to characterize their stability/unstability and to compute the periods of the equilibria.

The proposed theory is general enough, to be applied to a wide range of probes around planet or natural planetary satellites.

The BepiColombo mission is used to motivate our analysis and to provide specific numerical data to check our analytical results.

Finally, we also bring to the light that the coefficient J 2 is able to protect against the increasing of the eccentricity due to the Kozai-Lidov effect and the coefficient J 3 determines a shift of the equilibria.

Introduction

BepiColombo (MPO and MMO orbiters) is a joint European and Japanese space agencies mission aimed at studying the planet Mercury. The dynamics of the rotation of Mercury has been analyzed, for example, in D'Hoedt and [START_REF] D'hoedt | Planetary long periodic terms in Mercury's rotation: a two dimensional adiabatic approach[END_REF]; D 'Hoedt et al. (2010); [START_REF] Dufey | Latitudinal librations of Mercury with a fluid core[END_REF] . . . The MPO (Mercury Planetary Orbiter) will be brought into a polar elliptical orbit around Mercury with an inclination of 88-90 • , an eccentricity of 0.1632 and a semi-major axis of 3 394 km. The MMO (Mercury Magnetospheric Orbiter) will also be brought into a polar elliptical orbit with an eccentricity of 0.6679 and a semi-major axis of 8 552 km.

Actually polar orbits are very interesting for scientific missions to planetary satellites (with near polar low-altitude) or to planet (with high-eccentric high-altitude). The orbital dynamics of such space probes is governed by the oblateness (J 2 effect) of the central body around which the space probe is orbiting and the gravity field from the third body. A well-known effect of the third-body perturbation is the change in the stability of circular orbits related to orbit inclination. This effect is a natural consequence of the Kozai-Lidov resonance [START_REF] Kozai | Secular perturbations of asteroids with high inclination and eccentricity[END_REF][START_REF] Lidov | Evolution of the orbits of artificial satellites of planets as affected by gravitational perturbation from external bodies[END_REF]. Generally the final fate of such a satellite is the collision with the central body. However if the radius of the central body is small, large cyclic variations in eccentricity are possible without collision; for examples see [START_REF] Russell | Circulating eccentric orbits around planetary moons[END_REF]. Therefore the control of the orbital eccentricity leads to the control of the satellite lifetime. [START_REF] Scheeres | Stabillity analysis of planetary satellite orbiters: Application to the europa orbiter[END_REF] studied near-circular orbits in a model that included both the third body's gravity and J 2 . In addition [START_REF] San-Juan | Phase space structure around oblate planetary satellites[END_REF] studied orbit dynamics about oblate planetary satellites using a rigorous averaging method. [START_REF] Paskowitz | Design of science orbits about planetary satellites: Application to europa[END_REF] added the effect of the coefficient J 3 . These authors mainly focused their attention to an orbiter around planetary satellites especially for Europa. So they did not take into account the eccentricity of the third body and they detailed the nearcircular orbits. Recently [START_REF] Lara | Analytical theory for spacecraft motion about Mercury[END_REF] studied the planar elliptic restricted threebody problem including the J 2 and J 3 coefficients, using Deprit's perturbation method allowing to build a doubly averaged Hamiltonian and to provide the transformation between osculating and mean initial conditions. However they were not able to provide either the new stable horizontal equilibrium induced by the J 2 term, or the analytical formulation of the position of the equilibria and their periods.

Our purpose is to build an easy method to find a simplified Hamiltonian model, as simple as possible, which will reproduce the motion of probes orbiting an oblate central body also taking into account the third body effect with a non-null third body eccentricity. Especially we are looking for the conditions that give rise to frozen orbits. Frozen orbits are orbits that have orbital elements constant on average. These particular orbits are able to keep constant the eccentricity. Therefore in a neighbourhood of these orbits there is a stability area where even a limited control could be used to avoid the crash with the central body. Here we are interested in determining all generic frozen orbits, not only the circular equilibria. Notice that also periodic orbits can be used to find high fidelity orbits including the first order gravity effects and arbitrary order gravity field terms [START_REF] Lara | Computation of a Science Orbit about Europa[END_REF][START_REF] Russell | Long-lifetime Lunar Repeat Ground Track Orbits[END_REF].

Besides the oblateness (J 2 gravity term) of the central body and the gravity effect of the third body, our averaged model takes also into account the eccentricity of the orbit of the third body. Moreover let us observe that our results are given in closed form with respect to the eccentricity and inclination of the probe, namely we do not perform any power series expansion; therefore, our theory applies for arbitrary eccentricities and inclinations of the space probe, and is not limited to almost-circular orbits. We can thus conclude that the theory is general enough to be applied to a wide range of probes around a planet or around a natural planetary satellite and, can be formulated and presented in a general way that allows extension of the results to other cases.

The Mercury orbiter mission (BepiColombo) is used to motivate our analysis and to provide specific numerical data to check our analytical results.

We are able to provide the location of all frozen orbits and study their stability as a function of the involved parameters, using implicit equations and graphics. Finally we give the analytical expressions of the periods at the stable equilibria.

The analytical results are verified and confirmed using dedicated numerical simulations of the complete model.

To conclude, we discuss the effect of protection of J 2 on the increase of the eccentricity due to Kozai-Lidov effect and the apparition of an asymmetry caused by the addition of the coefficient J 3 .

These equilibria and their periods can be modified by including higher order terms of the gravity field. Except for resonant relations between the frequencies of the satellite and the central body, they will not change drastically the behavior of our phase space. Henceforth, the frozen orbits presented here provide good initial conditions to find periodic or frozen orbits of the full model.

Motivation: numerical exploration

For the purpose of our study, we consider the modeling of a space probe subjected to the influence of Mercury's gravity field (in the following sections Mercury will be denoted by "central body") and the gravitational perturbations of the Sun (noted "third body") as well as to the direct solar radiation pressure without shadowing effect. As a consequence the system of differential equations describing the probe motion is given by r = rpot + r⊙ + rrp ,

where rpot is the acceleration induced by the Mercury's gravity field which can be expressed as the gradient of the following potential

U pot (r, λ, φ) = - µ r + µ r ∞ n=2 n m=0 Rp r n P m n (sin φ)(Cnm cos mλ + Snm sin mλ) , (2)
where the quantities Cnm and Snm are the spherical harmonics coefficients of the hermeopotential. µ and Rp are respectively the gravitational constant and equatorial radius of Mercury. The quantities (r, λ, φ) are the hermeocentric spherical coordinates of the space probes. P m n are the associated Legendre polynomials. The acceleration r⊙ results from the gravity interaction with the Sun and can be expressed with respect to Mercury's centre of mass:

r⊙ = -µ ⊙ r -r ⊙ r -r ⊙ 3 + r ⊙ r ⊙ 3 , (3)
where r and r ⊙ are the hermeocentric coordinates of the space probe and of the Sun respectively. The quantity µ ⊙ is the gravitational constant of the Sun. This acceleration can also be expressed as the gradient of the following potential:

U ⊙ = µ ⊙ 1 r -r ⊙ - r • r ⊙ r ⊙ 3 . (4)
Regarding the direct solar radiation pressure, we assume an hypothetically spherical space probe with optical properties defined by a single scalar coefficient. The albedo of Mercury is ignored and Mercury's shadowing effects are not taken into account either. Similarly, the acceleration induced by the direct solar radiation pressure is modeled by

rrp = Cr Pr a ⊙ r -r ⊙ 2 A m r -r ⊙ r -r ⊙ , ( 5
)
where Cr is the adimensional reflectivity coefficient (fixed to 1 in the rest of the paper) which depends on the optical properties of the space object surface; Pr = 4.56 × 10 -6 N/m 2 is the radiation pressure for an object located at the distance of 1 AU from the Sun; a ⊙ is a constant parameter equal to the mean distance between the Sun and Mercury. The coefficient A/m is the area-to-mass ratio of the space object and will be fixed to 0.01m 2 /kg.

It is worth noting that we model the gravity potential of the central body only using the J 2 , C 22 and J 3 coefficients. In our implementation, we choose the highly accurate Solar System ephemeris given by the Jet Propulsion Laboratory (JPL) to provide the positions of the Sun and Mercury [START_REF] Standish | JPL planetary and lunar ephemeris[END_REF]. We adopt the variable step size Bulirsch-Stoer algorithm (see e.g. [START_REF] Stoer | Introduction to numerical analysis[END_REF] to numerically integrate the differential equation (1). Let us note that, for the purpose of validation, we also use a second numerical integrator DOP853 (an explicit Runge-Kutta method of order 8(5,3) with stepsize control due to Dormand & Prince [START_REF] Hairer | Solving ordinary differential equations I. Nonstiff problems[END_REF])).

In Figure 1 we report the results of a numerical integration of the system of equations (1) for a set of 19 600 orbits, propagated over a 200 year time span with a entry-level step size of 300 seconds. We consider a set of initial conditions defined by an eccentricity grid of 0.005 and a semi-major axis grid of 35 km, spanning the [2600, 7600] km range. The other fixed initial conditions are i 0 = 90 • for the inclination, Ω 0 = 67.7 • ω 0 = -2 • for the longitude of the ascending node and the argument of periherm, respectively; the mean anomaly at epoch is fixed to M 0 = 36.4 • on 14 September 2019. The area-to-mass ratio A/m = 0.01m 2 /kg. These values have been fixed by the initial conditions of BepiColombo mission found in [START_REF] Garcia | Bepicolombo mercury cornerstone consolidated report on mission analysis[END_REF].

Fig. 1 The eccentricity computed as a function of the initial eccentricity e 0 and the initial semi-major axis a 0 . The equations of motion include the central body attraction, the harmonics J 2 , C 22 , J 3 , the solar interaction as well as the perturbing effects of the solar radiation pressure (A/m = 0.01m 2 /kg). The eccentricity step is 0.005 and the semi-major axis step is 35 km. The initial conditions are i 0 = 90 We show the amplitude of the eccentricity (that is the difference between the maximum and minimum eccentricity reached during the integration) of each orbit in the left panel of Figure 1. For each orbit, using the Numerical Analysis of Fundamental Frequencies, for short NAFF [START_REF] Laskar | Secular evolution of the solar system over 10 million years[END_REF][START_REF] Laskar | Frequency map analysis and quasiperiodic decomposition[END_REF], we also calculate the fundamental frequency (noted ν) of the evolution of the eccentricity vector (k = e cos ω, h = e sin ω). We plot the logarithm of the second derivative (noted log | ∂ 2 ν ∂k ∂h |) of this frequency in the right panel of Figure 1, namely an indicator of the diffusion in the frequency space, hence also of the regularity of the orbit when it assumes small values. For more details concerning this use of frequency analysis, see [START_REF] Lemaître | A web of secondary resonances for large A/m geostationary debris[END_REF] where the frequency analysis has been used to study resonances in Geostationary Earth Orbits.

First, let us observe that the white zone in Figure 1 corresponds to orbits that crash onto central body's surface. Second we distinguish a curve where the variation of the eccentricity amplitude is null (dashed black line). On the right panel we plot log | ∂ 2 ν ∂k ∂h | (where ν is the fundamental frequency of the signal (k, h)) and we also distinguish on the left of the dashed black line (null-variation of eccentricity) a larger value of the log of the derivative that could correspond to a separatrix. These structures will be analyzed and explained using a simplified model, that takes into account the central body attraction with the J 2 harmonic coefficient and the third body gravitational effect. We observed that the solar radiation pressure (for further details on the effect of direct solar radiation onto BepiColombo, we refer to [START_REF] Garcia | Bepicolombo mercury cornerstone consolidated report on mission analysis[END_REF] and [START_REF] Lucchesi | The Non-Gravitational Perturbations impact on the BepiColombo Radio Science Experiment and the key rôle of the ISA accelerometer: direct solar radiation and albedo effects[END_REF]) and the C 22 coefficient do not play any role in these structures, hence these effects will be absent in the simplified model.

The Hamiltonian Formalism

The aim of this section is to introduce the Hamiltonian (11) already found in Tremaine et al (2009). Kepler's Hamiltonian describing the motion of a test particule orbiting an isolated point mass M is

H K = 1 2 v 2 - GM r = - GM 2 a ( 6
)
where G is the gravitational constant, r is the planetocentric position of the particule, v = ṙ, r = |r| and a is the semi-major axis of the particule.

To expand our potential, we use a vectorial method developed by Tremaine et al (2009) and used, for example, by Farago and Laskar (2010).

The quadrupole potential arising from an oblate planet ("central body") is (Eq. 2 with n = 2)

Φ J2 (r) = GM J 2 R 2 p r 3 P 2 (sin φ) = GM J 2 R 2 p r 3 P 2 cos π 2 -φ ¡ = GM J 2 R 2 p 2r 5 3(r •np) 2 -r 2
(7) where np is the unit vector oriented to the central body's spin axis (see Figure 2). M , Rp and J 2 are respectively the mass, the radius and the oblateness coefficient of the central body (planet or natural satellite).

We assume that r ≪ a 3b (where the subscript 3b is related to the third body), then the potential (4) can be formulated in terms of the Legendre polynomials (where Ψ is the planetocentric angle between the third body and the particule):

U 3b = GM 3b r 3b ∞ n≥2 r r 3b n Pn(cos Ψ ) (8)
Now, we only keep the second order Legendre polynomial (n = 2) and we average over the third body orbital period. So, we obtain the quadrupole in terms of the third body gravitational effect

Φ 3b (r) = GM 3b 4a 3 3b (1 -e 2 3b ) 3/2 3(r • n 3b ) 2 -r 2 (9)
where n 3b is the normal to the central body orbit. M 3b , a 3b and e 3b are respectively the mass, the semi-major axis and the eccentricity of the third body. This quadrupole term takes into account the eccentricity (e 3b ) of the third body (e.g. Sun or Jupiter) around the central body (planet or natural satellite). Let us stress the fact that [START_REF] Paskowitz | Orbit mechanics about planetary satellites[END_REF]; [START_REF] San-Juan | Phase space structure around oblate planetary satellites[END_REF]; [START_REF] Scheeres | Stabillity analysis of planetary satellite orbiters: Application to the europa orbiter[END_REF] do not include this eccentricity factor in their formulation, while for a Sun-Mercury-orbiter application, this will be an important contribution.

We then average over the Keplerian orbit of the test particule described by the following elements: a semi-major axis a, an eccentricity e, and an orientation specified by the unit vectors n along the angular momentum vector, u toward the pericenter and v = n × u. We have [START_REF] Brouwer | Methods of Celestial Mechanics[END_REF])

< r 2 > = a 2 1 + 3 2 e 2 , 1 r 3 = 1 a 3 (1 -e 2 ) 3/2 , < (r • u) 2 > = a 2 1 2 + 2e 2 , < (r • v) 2 > = a 2 1 2 -1 2 e 2 , (r • u) 2 r 5 = (r • v) 2 r 5 = 1 2a 3 (1 -e 2 ) 3/2 .
(10) where < > denotes the average over M , the mean anomaly of the orbit.

Let j ≡ Ô 1 -e 2 n , e = eu, τ = Ö GM a 3 t, ε J2 = J 2 R 2 p a 2 and ε 3b = M 3b a 3 M a 3 3b (1 -e 2 3b ) 3/2
where e is the eccentricity vector and ε J2 ≥ 0, ε 3b ≥ 0. We finally define a dimensionless (divided by GM/a) Hamiltonian

K ′ = - 1 2 + ε J2 4(1 -e 2 ) 5/2 1 -e 2 -3(j • np) 2 + 3ε 3b 8 5(e • n 3b ) 2 -(j • n 3b ) 2 -2e 2 , ( 11
)
That describes the secular equations of motion of a test particule around an oblate central body perturbed by the third body gravitational effect. Let us summarize the used assumptions:

1. the precession rate of the central body spin due to third body gravity is negligible; 2. the satellite is a massless test particule; 3. the third body is far enough from the central body such that the third body gravity can be approximated by a quadrupole; 4. the satellite is far enough from the central body so that the potential from the central body can be approximated as a monopole plus a quadrupole; 5. the perturbing forces (Φ J2 + Φ 3b ) are weak enough so that the secular equations of motion can be used to describe the orbital motion;

6. there are not resonant relations in mean motions between the frequencies of the satellite and the frequencies of the central body.

Let us remark that San-Juan et al. ( 2006) already studied the orbit dynamics about planetary satellites using an extensive averaging method based on the Lie transforms to obtain averaged equations involving higher orders whose result is the introduction of an asymmetry for direct and retrograde satellite. Our simplified model will not be able to capture this asymmetry because the resulting Hamiltonian (13) will be symmetric in the satellite inclination; thus direct and retrograde satellites will have the same behavior.

Let us now make some assumptions suitable in the case of a non-inclined central body orbit (e.g. Sun-Mercury-orbiter system or Jupiter-Europa-orbiter system). We hereby consider an equatorial third body, thus np = n 3b . We also set G = √ 1 -e 2 and H = G cos ı where j

• np = √ 1 -e 2 cos ı.
To eliminate an extra parameter, we divide the Hamiltonian by the coefficient ε J2 and we introduce the coefficient

γ γ = ε 3b ε J2 = M 3b M a 3 3b (1 -e 2 3b ) 3/2 a 5 J 2 R 2 p . ( 12
)
In Figure 2, we represent the geometry for the general problem (on the left) and for our simplified one (on the right). The averaged Hamiltonian is then

K ′ = ε J2 4G 3 1 -3 H 2 G 2 + 3ε 3b 8 5(1 -G 2 ) 1 - H 2 G 2 sin 2 ω -H 2 -2 + 2G 2 K ′ /ε J2 noted by K ⇐⇒ K = 1 4G 3 1 -3 H 2 G 2 + 3γ 8 5(1 -G 2 ) 1 - H 2 G 2 sin 2 ω -H 2 -2 + 2G 2 .
(13) This one degree of freedom Hamiltonian is independent of the ascending node Ω. If we take an inclined third body, this ascending node takes place in the Hamiltonian and we obtain a two degree of freedom problem which is non integrable. If we take γ = 0 (ε 3b = 0, namely we only take into account the oblateness effect), we have the well-known circular dynamics of the eccentricity vector due to the J 2 coefficient with an elliptical fixed point in the semi-equinoctial elements (k

, h) = ( √ 1 -G 2 cos ω, √ 1 -G 2 sin ω).
If we take γ → ∞ (ε J2 = 0 i.e. only the third body contribution does matter), we find the Kozai-Lidov Hamiltonian which can be found in a similar formalism in [START_REF] Paskowitz | Orbit mechanics about planetary satellites[END_REF] (with e 3b = 0). The Hamiltonian (13) (with e 3b = 0) can also be found in [START_REF] Scheeres | Stabillity analysis of planetary satellite orbiters: Application to the europa orbiter[END_REF].

Table 1 Connection between the semi-major axis (km) of the probe around the central body and the coefficient γ. The rows "Min.", "Missions" and "Hill" give the values of γ with respect to the semi-major axis respectively equal to the radius of the central body, to one space mission and to the radius of Hill's sphere. In the last columns, we take two particular values of γ (1/7 and 0.5) that will play a relevant role in the next sections. For Europa the minimum value of γ is 1.153, so this cell is empty. For illustration, we show in Table 1 and draw in Figure 3 the value of the coefficient γ with respect to the semi-major axis for a probe around a terrestrial planet and around Europa. This coefficient can be related to other parameters used in the literature. For example, it can be linked to the coefficient β in [START_REF] San-Juan | Phase space structure around oblate planetary satellites[END_REF], to the coefficient ǫ used in [START_REF] Scheeres | Stabillity analysis of planetary satellite orbiters: Application to the europa orbiter[END_REF] or to the coefficient J2 in [START_REF] Lara | Analytical theory for spacecraft motion about Mercury[END_REF].

Min

Secular Equations of Motion

From the doubly averaged Hamiltonian (13), we obtain the equations of motion:

Ω = -H ´3 2G 5 + 3γ 8 10 G 2 (1 -G 2 ) sin 2 ω + 2 µ (14) Ḣ = 0 (15) ω = 3γ 4 5 H 2 G 3 -G sin 2 ω + 2G + 3 4G 4 5 H 2 G 2 -1 (16) Ġ = - 15γ 4 (1 -G 2 ) 1 - H 2 G 2 sin ω cos ω . ( 17
)
Developing these equations in eccentricity up to second order, we can obtain the equations of [START_REF] Scheeres | Stabillity analysis of planetary satellite orbiters: Application to the europa orbiter[END_REF]. In the following, we will adopt a complementary approach, keeping functions of eccentricity and inclination, without power series developments, in such a way that our results hold for arbitrary eccentricities and inclinations.

From the previous set of equations, we observe that H 2 = G 2 cos 2 ı is a constant of motion as in the Kozai-Lidov effect [START_REF] Kozai | Secular perturbations of asteroids with high inclination and eccentricity[END_REF][START_REF] Lidov | Evolution of the orbits of artificial satellites of planets as affected by gravitational perturbation from external bodies[END_REF]. Besides, let us remark that 0 ≤ G ≤ 1, thus 0 ≤ H ≤ G ≤ 1 and moreover γ > 0. The first equation ( 14) is equal to zero only for H = 0 corresponding to exact polar inclination. Moreover the ascending node does not affect any of the other orbital elements. The last equation (17) equals to zero for G = 1, ω = kπ/2, k ∈ N or H = G, namely ı = 0 • , that is the planar case. The third equation ( 16) could equal to zero for ω = 0, π or ω = ±π/2. We analyze these equations in next section to find the equilibria.

Frozen Orbit Solutions

A frozen orbit is characterized by no secular change in orbital eccentricity and argument of pericenter. It has constant values of e, ı and ω on average, this results in fixed geometrical size and locations, apart from short period oscillations.

We have already observed that equilibria appear when G = 1 or ω = 0, π or ω = ±π/2. We separatly deal with these three different cases. For each of them, we give the number of equilibria, the conditions of existence and we calculate their stability.

We do not deal with the singularity G = 0 (⇔ e = 1) because it corresponds to an escape of the orbiter. We will show that the equilibrium G = 1 always exists. So to begin, we deal with the non-circular case G = 1 (eccentricity = 0).

5.1 Non-circular case G = 1 (eccentricity = 0)

5.1.1 Vertical equilibria -Kozai-Lidov equilibria: cos ω = 0 ⇔ ω = ±π/2
The conditions to simultaneously equal to zero the equations ( 16) and ( 17) is:

H 2 = G 2 5 1 + 3G 5 γ 1 + G 3 γ cos ω = 0 (18)
Because 0 ≤ G < 1 then the first equation of ( 18) implies that

H 2 < 1 + 3γ 5γ + 5 . ( 19
)
Let us observe that this is also the value of H2 for which one real root does exist. If this condition is violated then no real root exists.

Actually we determine1 a region 2 given by the implicit equation 864 000

H
)
where it is possible to find three real roots. We will show that these three reals roots appear for eccentricities larger than 0.996 59. Being a case close to an escape of the orbiter, we will leave to section 7.3.3 a discussion of this "local deformation".

If the oblateness term is neglected (ε J2 ≃ 0 ⇔ γ → ∞), the existence condition becomes independent of the physical parameter and reduces to sin 2 ı < 2 5 or arccos

Õ 3 5 ≃ 39.23 • ≤ ı ≤ 144.77 • which corresponds to Kozai-Lidov critical inclina- tion.
We also analyze the stability of these equilibria (18). The Jacobian of the Hamiltonian (13) evaluated at the equilibrium (18) (noted by | Eq.( 18) or G kl being the value of G at the Kozai-Lidov equilibrium) is given by:

∂ 2 K ∂G 2 ¬ ¬ ¬ ¬ Eq.(18) = 3 2G 5 2 -15 H 2 G 2 - 9γ 4 1 + 5 H 2 G 4 ¬ ¬ ¬ ¬ Eq.(18) = 3 4G 5 kl 1 1 + γG 3 kl -2 + γG 3 kl -21γG 5 kl -12γ 2 G 8 kl ∂ 2 K ∂ω 2 ¬ ¬ ¬ ¬ Eq.(18) = - 15 4 γ(1 -G 2 ) 1 - H 2 G 2 ¬ ¬ ¬ ¬ Eq.(18) = - 3 2 γ(1 -G 2 kl ) 2 -G 5 kl γ 1 + G 3 kl γ ∂ 2 K ∂G ∂ω ¬ ¬ ¬ ¬ Eq.(18) = ∂ 2 K ∂ω ∂G ¬ ¬ ¬ ¬ Eq.(18) = 0. (21)
In the equations ( 21), the term

∂ 2 K ∂ω 2 ¬ ¬ ¬
Eq.( 18)

is always strictly negative (if G < 1). Then the equilibrium is a stable point if

1 G 5 2 -15 H 2 G 2 - 3γ 2 1 + 5 H 2 G 4 ¬ ¬ ¬ ¬
Eq.( 18)

< 0 ⇐⇒ -2 -21γG 5 kl + 15γG kl H 2 < 0 .

(22) This equation ( 22) is always satisfied for all γ > 0, H 2 < 1+3γ 5γ+5 and G kl < 1 (e kl > 0). Therefore, for these conditions, we have two opposite stable points at ω = ±π/2 and G such that

H 2 = G 2 5 1+3G 5 γ 1+G 3 γ .
For the inclination ı = 90 • (H 2 = 0) the equilibrium exists for the particular value of G = 0 (e = 1). This case is only theoretical and should not be considered, because it would correspond to an escape of the orbiter. In Figure 4 we give the location of the equilibria G (18) in the parameter space (γ,H 2 ). The conditions to simultaneously equal to zero the equations ( 16) and ( 17) are:

H 2 = G 2 5 (1 -2G 5 γ) sin ω = 0. (23)
Using "Le théorème d'algèbre de Sturm" (Sturm 1835) we calculate the number of roots (G) in the range 0 ≤ G < 1 of the equation ( 23) as a function of the parameters γ and H 2 . For γ > 0 this equation has

• one real root, equal to 0 if H 2 = 0 and γ < 1/2.

• three real roots (one equal to 0 and the other two opposite) if H 2 = 0 and γ ≥ 1/2.

• three real roots (one equal to 1 and the other two opposite) if 0 < H 2 < 1-2γ 5 ; • five real roots (one equal to 1 and the other ones opposite two by two) if γ ≥ 1/7 and 1-2γ 5 < H 2 < (7γ) -2/5 7 ; • one real root equal to 1 otherwise.

In Figure 5, we give the location of the equilibria G (23) in the space (γ,H 2 ). The particular case G = 1 will be treated in the next section. We can also analyze the stability of these equilibria (23). The Jacobian of the Hamiltonian (13) evaluated at the equilibrium (23) (noted by | Eq.( 23) or G hor , being G hor the value of G at the equilibrium) is given by:

∂ 2 K ∂G 2 ¬ ¬ ¬ ¬ Eq.(23) = 3 2G 5 2 -15 H 2 G 2 + 3γ 2 ¬ ¬ ¬ ¬ Eq.(23) = 3 2G 5 hor -1 + 7γG 5 hor ∂ 2 K ∂ω 2 ¬ ¬ ¬ ¬ Eq.(23) = 15 4 γ(1 -G 2 ) 1 - H 2 G 2 ¬ ¬ ¬ ¬ Eq.(23) = 3 2 γ(1 -G 2 hor )(2 + G 5 hor γ) ∂ 2 K ∂G ∂ω ¬ ¬ ¬ ¬ Eq.(23) = ∂ 2 K ∂ω ∂G ¬ ¬ ¬ ¬ Eq.(23) = 0. ( 24
)
In the equations ( 24), the term

∂ 2 K ∂ω 2 ¬ ¬ ¬
Eq.( 23)

is always strictly positive (if G < 1). Then the equilibrium is a stable point if γ > 1 G 5 15 H 2 G 2 -2 ⇐⇒ G 5 hor > 1 7γ . (25)
Using this equation at the equilibrium (23), we obtain conditions for stability of the stable point (G = 1 ⇔ e = 0)

1 -2γ 5 < H 2 < 1 7 1 7γ 2/5 and 1 7 ≤ γ . (26)
So the condition to have an unstable equilibrium is given by

γ < 1 7 or H 2 < 1 -2γ 5 or H 2 > 1 7 1 7γ 2/5 . ( 27
)
In the Figure 5, we notice that when both unstable and stable equilibrium exist, the unstable equilibrium always appears for a value of G lower than that of the stable point (i.e. for a value of e greater than the one for the stable point). For the case G = 1, we can use a canonical transformation to cartesian coordinates

x = Ô 1 -G 2 sin ω y = Ô 1 -G 2 cos ω (28)
The new Hamiltonian is therefore

K = 1 4 1 (1 -x 2 -y 2 ) 3/2 - 3H 2 (1 -x 2 -y 2 ) 5/2 (29) + 3γ 8 5x 2 1 - H 2 1 -x 2 -y 2 -H 2 -2x 2 -2y 2
for which it is obvious that (0, 0) is always an equilibrium point whose stability can be studied computing the second derivatives and evaluate them at this equilibrium:

∂ 2 K ∂x 2 ¬ ¬ ¬ x=0=y = 3 4 (1 -5H 2 ) + 3γ 4 (3 -5H 2 ) ∂ 2 K ∂y 2 ¬ ¬ ¬ x=0=y = 3 4 (1 -5H 2 ) - 3γ 2 ∂ 2 K ∂x ∂y ¬ ¬ ¬ x=0=y = ∂ 2 K ∂y ∂x ¬ ¬ ¬ x=0=y = 0 . (30)
So, the condition to have a stability point at x = 0 = y is

H 2 < 1 -2γ 5 or H 2 > 1 + 3γ 5γ + 5 ; ( 31
)
and thus the condition to have an unstable point at

x = 0 = y is 1 -2γ 5 < H 2 < 1 + 3γ 5γ + 5 . ( 32
)
5.3 Summary of the phase space

In this section we summarize the topology of the various possible phase spaces as a function of the parameters. We draw (Fig. 6) the bifurcation lines (conditions 19, 26 and 31) in the parameter space (γ,H 2 ). This bifurcation diagram is equivalent to the upper part of the bifurcation diagram in San-Juan et al. ( 2006) but here we draw the bifurcation lines in the general (not linked to a particular central body) space (γ,H 2 ). Let us observe that the bifurcation diagram presented in figure 3 of [START_REF] Lara | Analytical theory for spacecraft motion about Mercury[END_REF] shows only the conditions (31) because authors were only interested in circular orbits (when J 3 = 0). The H 2 = (7γ) -2/5 /7 line stops at the limit γ = 1/7. For this value, this curve coincides with the H 2 = (1 -2γ)/5 condition. For the Jupiter-Europa-orbiter system, the minimum value of γ is 1.153 (Tab. 1). Therefore the phase spaces (A) and (E') do not exist.

Fig. 6 Bifurcation lines and regions in the parameter space (γ,H 2 ). In the regions (B) and (C) we have the same number and the same stability of the equilibria but the phase space is topologically different. These regions are separated by the dashed line implicitily given by equation ( 33). For each region there are stable and/or unstable point(s). To visualize these points, see Fig. 7.

The region (E') and (F) in magenta color correspond to exact polar orbits (ı = 90 • thus H = 0).

For each region, we attribute a letter and we draw (Fig. 7) a generic contour plot of the Hamiltonian (13) in the (k, h, ı) physical space. We recall that the motion of the inclination ı is given by the conservation of the first integral H = G cos ı. We also draw the projection of these phase spaces in the semi-equinoctial elements space (k, h) = ( √ 1 -G 2 cos ω, √ 1 -G 2 sin ω). In this phase space, it is easier to bring to the fore the stable (green point) and unstable (red cross) equilibria. The (E') phase space is trivial, containing only concentric circles in the ı = 90 • plane, so we do not repoduce it.

In the Figure ( 7), we notice that the maximum inclination is always reached at e = 0. This is explained by the relation H 2 = √ 1 -e 2 cos ı. This last relation also gives a maximum bound onto the eccentricity: e ≤ √ 1 -H 2 . Therefore there are some values of H for which the phase space is visibly restricted in eccentricity. Beyond this eccentricity, the motion is physically impossible.

k, h) = ( √ 1 -G 2 cos ω, √ 1 -G 2 sin ω).
The green point and red cross are respectively the stable and unstable points. In the polar projection, the radius is the eccentricity e and the angle is the pericenter ω in degrees.

Let us observe that the region near the stable equilibria allows designers to control the variation of the eccentricity even for high eccentricity. We also remark that there are "dangerous" portions of phase space such as the region around the γ = H 2 = 1/7 or near of the (B)-(C) transition. In these regions the dynamics (in a full model) could strongly change for a small variation of (H 2 , γ) or (e, ω).

The transition between (B) and (C) phase spaces arises when the energy of the separatrix at the (0, 0) equilibrium is equal to the energy of the unstable exterior horizontal equilibrium. This condition gives a new "fictitious" bifurcation line (dashed line in Figure 6) in the parameter space (γ, H 2 ). To find this line, we evaluate the Hamiltonian (13) at the unstable equilibrium H 2 = G 2 5 (1-2G 5 γ) (Eq.23 and condition 27) and we denote this value by K 1 . Afterward, we evaluate the Hamiltonian (29) at the unstable equilibrium (0, 0) (condition 32) and denote the result by K 2 . We now assume these two equilibria have the same value of Hamiltonian K and of H 2 . Then we can replace H 2 by G 2 5 (1 -2G 5 γ) in K 2 and we impose the equality between K 1 and K 2 . Therefore we obtain the condition

γ = 2 + 5G 3 hor + 3G 5 hor 6G 10 hor + 15G 3 hor -21G 5 hor where H 2 = G 2 hor 5 (1 -2G 5 hor γ) ( 33
)
where G hor is the unstable horizontal equilibrium i.e. H 2 = G 2 hor 5

(1 -2G 5 hor γ). We plot this implicit condition (33) in Figure 6 with a dashed black line. This line joins the "(1 -2γ)/5" and "(7γ) -2/5 /7" lines at the (γ = 1/7, H 2 = 1/7) point.

For the particular case γ → 0 (J 2 effect only), we obtain, for all H 2 , a phase space with circular motion of the eccentricity. We see that near to value H 2 = 1/2 (corresponding to the Molniya3 critical inclination equal to 63.43 • with G = 1; for more informations about critical inclination for a massive satellite refer to [START_REF] Breiter | Critical inclination in the main problem of a massive satellite[END_REF], the phase spaces (A), (D) and (E) always exist until γ becomes exactly equal to 0. In the opposite case, γ → ∞ (third body effect only), the curve H 2 = (7γ) -2/5 /7 converges to 0 and the curve H 2 = (1 + 3γ)/(5γ + 5) converges to 3/5 (corresponding to the Kozai-Lidov critical inclination equal to 39.23 • with G = 1). Then only the following phase spaces are realizable: (E) (for H 2 > 3/5), (D) (for 0 < H 2 < 3/5) and (F) (for H 2 = 0) with (F) that degenerates to an unstable point at the center. These three phase spaces will be shown in Figure 12.

In Figure 8, we show how the stable and unstable equilibria evolve, appear and disappear in each region and during the transition between the regions. We take a vertical section in Figure 6 at γ = 0.4. This section crosses the regions (A), (C), (B), (D) and (E). We draw the value of the eccentricity at the stable (solid color lines) and unstable (dashed color lines) equilibria with respect to H 2 . The vertical dashed black lines mark the boundary of the regions. The numbers give the number of equilibria with this value of e. For example, 2 in magenta dashed line means that there are two unstable equilibria with the same value of e, respectively for ω = 0 and ω = π.

At the transition between (A) and (C), the central (e = 0) stable point bifurcates in two horizontal stable points (e = 0 and ω = 0, π) and one unstable point (e = 0). At the transition between (B) and (D), the two unstable and the two stable horizontal (ω = 0, π) equilibria converge to the same value of e and cancel out. At the transition

A) (C) (D) (E) 2 2 2 2 2 1 1 1 1 2 (B)
Fig. 8 For a vertical section (γ = 0.4) in the Figure 6, value of the eccentricity at unstable (dashed color lines) and stable (solid color lines) equilibria with respect to H 2 . The numbers give the number of equilibria for each curve.

between (D) and (E), the two stable Kozai-Lidov equilibria (ω = ±π/2) come close to 0 and cancel out with the central unstable equilibrium to give one central stable equilibrium. We remark that the transition between (C) and (B) is not characterized by a change of the equilibria.

Period at the equilibrium

We are now interested in the period of the eccentricity vector at the equilibrium. This will be done by linearizing in a neighborhood of the equilibrium. Then the Hamiltonian close to the equilibrium is given by (the subscript eq. means "evaluated at the equilibrium"):

K = Keq. + ∂K ∂G ¬ ¬ ¬ eq. ßÞ =0 (G -Geq.) + ∂K ∂ω ¬ ¬ ¬ eq.
ßÞ =0

(ω -ωeq.)

+ 1 2 ∂ 2 K ∂G 2 ¬ ¬ ¬ eq. (34)
ßÞ not.

= a (G -Geq.

ßÞ not. = X ) 2 + ∂ 2 K ∂G ∂ω ¬ ¬ ¬ eq. ßÞ =0 (G -Geq.)(ω -ωeq.) + 1 2 ∂ 2 K ∂ω 2 ¬ ¬ ¬ eq.
ßÞ not.

= b

(ω -ωeq.

ßÞ not. = Y ) 2 K = Keq. + aX 2 + bY 2 . ( 35
)
This is an harmonic oscillator that can be expressed in action-angle variables (ψ, J) defined as (at a stable equilibrium, we have ab > 0):

X = 4 Ö b a √ 2J cos ψ and Y = 4 Ö a b √ 2J sin ψ . ( 36
)
Then the frequency at the equilibrium is given by

ψ = ∂K ∂J = 2 √ ab = × ∂ 2 K ∂G 2 ¬ ¬ ¬ eq. ∂ 2 K ∂ω 2 ¬ ¬ ¬ eq. . ( 37
)
Using the equation ( 37), the periods (τ ) at the stable equilibria are given by:

for horizontal equilibria: G such as stability (26

) J2 G 5 2 -15 H 2 G 2 + ε 3b ε 3b (1 -G 2 ) 1 -H 2 G 2 ; ( 38
)
for vertical (Kozai-Lidov) equilibria: G such as Equation ( 18) and condition of stability ( 19)

τ = 4π 3 Ö -5 2 ε J2 G 5 2 -15 H 2 G 2 - 3ε 3b 2 1 + 5 H 2 G 4 ε 3b (1 -G 2 ) 1 -H 2 G 2 ; ( 39
)
for central equilibrium (e = 0): G = 1 with condition of stability ( 31)

τ = 8π 3 Ö ε J2 (1 -5H 2 ) + ε 3b (3 -5H 2 ) ε J2 (1 -5H 2 ) -2ε 3b . ( 40
)
We remind that γ = ε 3b /ε J2 and that the equations are dimensionless. Then the periods at the equilibria are given by Teq. = Õ a 3 GM τeq.. For example, we apply these formula to a Mercury orbiter. The values for Mercury are a 3b = 57 909 176.0 km, e 3b = 0.205 630 69, J 2 = 6.0 × 10 -5 [START_REF] Anderson | The mass, gravity field, and ephemeris of mercury[END_REF] and Rp = 2 439.7 km [START_REF] Seidelmann | Report of the IAU/IAG Working Group on cartographic coordinates and rotational elements: 2006[END_REF]. In Figure 9 we plot the periods at the equilibria respectively for the three cases:

on the left panel, the periods at the stable equilibrium with respect to the values of γ and H 2 . The color code indicates the period of the fundamental frequency at the equilibrium; -on the right panel, the location of the stable equilibrium in the phase space (a, e, ı) with the period in the color scale.

The color code is the same for the left and right panels and it is truncated at the value of 100 years. For a larger period, we use the black color. We notice that the time of space missions is very short compared with those periods. However these long periods are interesting because they allow to keep constant eccentricity during the mission time.

6 Comparison of analytical and numerical solutions

Comparison for all inclinations

The analytical results of the simplified model described above are checked using a precise numerical integration of the complete set of equations of motion (1). For our test, we use Mercury's orbiter mission profile, which nominally puts the spacecraft into a high eccentric polar orbit. Numerical integrations were performed with the Bulirsch-Stoer [START_REF] Stoer | Introduction to numerical analysis[END_REF] integrator. We reproduce hereby a few characteristic plots of the numerical simulations to confirm our analytical theory (see Figure 10). Similar results have been obtained for a wide range of initial frozen orbit conditions.

Figure 10 shows a very good agreement between analytical results and numerical simulations over a large time span (180 years). ) and (0, 0) equilibrium (with conditions 31) respectively in the upper, center and lower panels. The color code indicates the period (truncated to 100 years) of the fundamental frequency at the equilibrium. On the left panels, the period with respect to the parameters (γ, H 2 ). On the right panels, the location of the stable equilibrium in the physical space (a, e, ı) with its period. For the equilibrium (0, 0), e is always equal to 0. The analytical model is based on Equations (14, 16, 17). We plot the numerical integrations with continued lines and the analytical results with dashed lines. The time of the integrations is equal to 180 years. In the right panel, the numerical integration leads to a crash onto the planet after 40 years.

Comparison for polar inclination and explanation of the preliminary numerical results

In the Figure 11, we present a graphical comparison between numerical integration and analytical results (contour plots of the Hamiltonian ( 13)) for an exact polar inclination. We see that the analytical theory is very close to the numerical integration for all initial eccentricities. We also notice that the addition of C 22 modify much the motion. Then the initial conditions of frozen orbits found by our method give good initial conditions to find frozen orbits in a full model.

In the right panel, we show two solutions close to the libration point and we see that, the closer the motion is to the libration equilibrium, the more the numerical integrations show a discrepancy with respect to the analytical results for the periherm libration: the frozen orbit of the analytical model shows no changes in eccentricity and argument of pericenter. On the contrary, the numerical orbit has short period oscillations but constant mean values of e and ω. For indication, the analytical and numerically calculated periods at this equilibrium are given in Table 2.

Figure 11 allows us to explain the behaviors already seen in our preliminary numerical exploration (Fig. 1). In Figure 11, we take different orbits with a semi-major axis equal to 6 000 km corresponding to a vertical section in Figure 1. Then, on this section, we take some values of the eccentricity such that:

for e near to 0, in Fig. 1, we see a large value of the amplitude of variation of the eccentricity approximatively equal to 0.5 and a high value of the second derivative. In Fig. 11, for e equal to 0, we are on the separatrix. Therefore the eccentricity increases (roughly until 0.5) and a little shift of the initial eccentricity causes a high difference of the frequency. Thus the second derivative of the frequency is large; -for e close to 0.37, in Fig. 1, we see that the amplitude of variation of the eccentricity decreases until 0. In Fig. 11, at e = 0.37, we find the stable point where the eccentricity is equal to a constant; -when e moves away from 0.37 to 0.5, in Fig. 1, we see that the amplitude of variation of the eccentricity increases from 0 to 0.5 and for e = 0.5, the amplitude of variation of the eccentricity is maximal and the value of the second derivative is large. In Fig. 11, moving away from the equilibrium (e = 0.37) toward the separatrix (e ≃ 0.5) we encounter larger and larger variations in e; -for e near to 0.55, in Fig. 1, we see that the amplitude of variation of the eccentricity is smaller than for e ≃ 0.5. In Fig. 11, for e ≃ 0.55, the pericenter circulates and the maximum of the amplitude of variation of the eccentricity is roughly equal to 0.55 -0.11 = 0.44. -in Fig. 1, moving along the line e = 0, we pass from the region (F) to the region (E') at 5 577 km (Tab.1). For semi-major axis smaller than a = 5 577 km, we do not cross any separatrix and the amplitude of variation of the eccentricity is small.

Frequency comparison

To obtain a second independent validation of our analytical model, we numerically compute, using the NAFF algorithm [START_REF] Laskar | Secular evolution of the solar system over 10 million years[END_REF][START_REF] Laskar | Frequency map analysis and quasiperiodic decomposition[END_REF], the period of the numerical solutions of the full system (1) obtained through numerical integration, and we compare it with the period of the equilibrium points of the simplified model.

Table 2 provides a summary of these comparisons. We can observe a very good agreement between the two methods. Some small differences can be explained as follows:

the exact equilibrium in the doubly averaged system is not the exact equilibrium in the full numerical model; the full numerical model contains short period terms which disturb the long period dynamics.

The closeness of the periods found analytically as well as numerically shows also that our method can provide a very good approximation of periods and location of frozen orbits even in the full osculating model.

Discussions

7.1 J 2 : the protector The aim of this section is to describe the protection mechanism of the coefficient J 2 on the increase of the eccentricity. We recall that our Hamiltonian (13), once we set the coefficient ε J2 = 0, reduces to the Kozai-Lidov Hamiltonian:

K kl = 3 ε 3b 8 5(1 -G 2 ) 1 - H 2 G 2 sin 2 ω -H 2 -2 + 2G 2 . ( 41
)
In the upper panels of Figure 12, we draw the possible phase spaces of this Hamiltonian.

For comparison, we plot in the lower panels analogous phase spaces with Mercury's J 2 = 6 × 10 -5 . In the right panels (H 2 > 3/5) we have a similar behavior of our (E) case (Fig. 7). For the exact polar orbits (H 2 = 0 in the left panel of the Fig. 12), in the Kozai-Lidov Hamiltonian (top left panel in Fig. 12), all the probes are ejected: the eccentricity always grows up to 1. Instead, with the addition of the coefficient J 2 we have (Fig. 7) the phase space (E') or (F) (in bottom left of Fig. 12 we have the E') where it is possible that the eccentricity does not increase or that it remains at a fixed value. In the middle case (0 < H 2 < 3/5) we see that for an initial pericenter close to 0, the eccentricity increases (top center panel in Fig. 12). Instead, in our case (J 2 = 0), the phase spaces (A), (B), (C) and (E) (Fig. 7) show that it is possible to find initial condition (other than ω ≃ ±π/2) where the increasing of the eccentricity is naturally controlled. Nevertheless, with the addition of the coefficient J 2 , we can keep the vertical eight-shape (phase space (D) in Fig. 7 or in bottom center in Fig. 12) but with a lower increase of eccentricity.

The J 2 acts as a protection mechanism against the increase of the eccentricity due Kozai-Lidov effect. This mechanism also appears for planets in tight binary systems [START_REF] Saleh | The stability and dynamics of planets in tight binary systems[END_REF], where the general relativistic effects become dominant and can cause the periastron to precess on very short timescales. Therefore this precession can lead to the suppression of Kozai oscillations.

Fig. 12 In the upper panels, all possible phase spaces for with respect to the values of H 2 . For comparison, in the lower panels analogous phase spaces with J 2 = 6 × 10 -5 adapted to the case of Mercury.

Local deformation of the Kozai-Lidov equilibrium

We have seen that the condition to get the Kozai-Lidov equilibrium is (Eq. 19)

H 2 < 1 + 3γ 5γ + 5 . ( 42
)
Actually there is a region, where it is possible to find three real roots for G, as a function of H 2 and γ. The conditions to have these three real roots are given by: (43) We draw the solutions of this equation, denoted by KL 3 , that demarcates the region denoted (G), on the left panel of the Figure 13. Let us observe that this condition is verified for large values of γ (γ ≥ 7203 √ 3/2) and for very small values 2 (H 2 ≤ 1/3087). An example of the phase space is plot in Figure 13 in the middle panels. In this region, the vertical Kozai-Lidov stable equilibrium bifurcates in two stable and one unstable vertical Kozai-Lidov equilibria producing thus a local deformation of the Kozai-Lidov equilibrium. We show an example of these three equilibria in the right panels of the Figure 13. Initial conditions close to these equilibria (external orbit in the right panels of Fig. 13) give rise to orbit librating around this set of three equilibria.

KL 3 ≡ 864
It is possible to find that this bifurcation appears, in the (G) region, for a value of G smaller than √ 3/21 ≃ 0.082 478 6 corresponding to a value of the eccentricity e larger than √ 438/21 ≃ 0.996 59. Recalling the formula H = G cos ı, we obtain a minimal inclination of 87.27 • . Nevertheless for the Sun-Mercury-orbiter application these conditions lead orbiters to crash onto Mercury.

J 3 discussion

In [START_REF] Paskowitz | Design of science orbits about planetary satellites: Application to europa[END_REF] and [START_REF] Lara | Analytical theory for spacecraft motion about Mercury[END_REF] the authors included the J 3 (the "pear shape") of the central body in their system. [START_REF] Paskowitz | Design of science orbits about planetary satellites: Application to europa[END_REF] noticed that the coefficient J 3 caused an asymmetry between the solutions of the frozen orbits for ω = ±π/2 but they did not explain the reasons of this beavior. [START_REF] Lara | Analytical theory for spacecraft motion about Mercury[END_REF] gived the values of the pericenter at the "horizontal" equilibria and a region (e, H) where frozen orbits may exist with ω = ±π/2.

The potential arising from a central body with a J 3 = 0 is given by

Φ J3 (r) = GM J 3 R 3 p 2r 7 (r • np) 5(r • np) 2 -3r 2 . ( 44
)
The averaged Hamiltonian

3 GM J 3 R 3 p 2 a 4 (1 -e 2 ) 5/2 e sin ω sin ı 1 - 5 4 sin 2 ı . (45)
It is easy to understand that the north-south asymmetrical form of J 3 induces an additional term proportional to sin(ı) (which is not an even function) in the Hamiltonian (45). The odd powers of (r • np) are responsible for emergence (after average) of the (sin ω sin ı) coefficient.

Using our variables G = √ 1 -e 2 , H = G cos ı, we can define the dimensionless (divided by GM/a) potential that we can add to the Hamiltonian (13):

J 3 R 3 p a 3 ßÞ not. = ε J3 3 8G 8 Ô 1 -G 2 sin ω Ô G 2 -H 2 (5H 2 -G 2 ). ( 46
) Introducing the coefficient δ = ε J3 ε J2 = J 3 Rp J 2 a
, the equations of motion (17 and 16) can be rewritten in compact form as follows:

Ġ = F 1 (G, H, γ) sin ω cos ω + F 2 (G, H, δ) cos ω ω = F 3 (G, H, γ) + F 4 (G, H, γ) sin 2 ω + F 5 (G, H, δ) sin ω (47)
where the functions F 1 , F 3 and F 4 can be easily identified in equations ( 16) and ( 17). The functions F 2 and F 5 come from the J 3 effect and they are proportional to δ. Let us observe that the addition of J 3 causes an asymmetry in the frozen orbit solutions not present before. Indeed, for ω = π/2 the condition of equilibrium is given by

F 3 (G, H, γ) + F 4 (G, H, γ) + F 5 (G, H, δ) = 0 (48)
whereas for ω = -π/2 the condition of equilibrium is given by

F 3 (G, H, γ) + F 4 (G, H, γ) -F 5 (G, H, δ) = 0. (49)
Then, for a small coefficient δ, the asymmetry is not important. However for a large value of this coefficient, the asymmetry could be important until the elimination of one of two equilibria.

7.3.2 Horizontal equilibria.

For horizontal equilibria, the condition of equilibrium ( Ġ = 0) becomes:

F 1 (G, H, γ) sin ω + F 2 (G, H, δ) = 0 ⇐⇒ sin ω = -F 2 /F 1 not. = -ǫ. (50)
Then the "horizontal" equilibria appear for non-zero values of the pericenter ω = -ǫ and ω = π + ǫ. The condition to obtain ω = F 3 (G, H, γ) + F 4 (G, H, γ)ǫ 2 -F 5 (G, H, δ)ǫ = 0, (51) that induces a shift in the equilibrium in G and ω variables with respect to the case "J 2 + third body". It is possible, in a similar way as before for J 2 , to calculate the new values of the pericenter at the equilibria by solving the first equation ( 47). Then Ġ vanishes for ω = ±π/2 (see previous section) or for

ω = arcsin - δ γ 5H 2 -G 2 10G 6 √ 1 -G 2 √ G 2 -H 2 (52)
with an existence condition on ω ∈ R bounded by

H 2 = G 2 5 1 -2G 5 γ 5G 5 γ(1 -G 2 ) ± 2 Õ 5 (1 -G 2 ) (5G 10 γ 2 (1 -G 2 ) + 4δ 2
) .

(53) These solutions are equivalent to the equation ( 44) in [START_REF] Lara | Analytical theory for spacecraft motion about Mercury[END_REF]. If we insert this solution in the second equation ( 47), ω vanishes for

H 2 = - 10G 7 γ -6δ 2 ± Ô 100G 14 γ 2 + 20G 7 γδ 2 + δ 4 -280γ 2 G 12 δ 2 G 2 35δ 2 . ( 54
)
Consequently with the method described in section 5.4 we can as well calculate the stability and the period of these equilibria.

Modifications of the phase space

For illustration, in Figure 14, we draw the contour plots of the new Hamiltonian for different values of J 3 (or for different values of δ). We see that when the δ coefficient increases (in absolute value), the vertical equilibrium goes down while the horizontal equilibrium goes below the "line sin ω = 0". We point out that, from some values of δ, the equilibrium ω = -π/2 disappears (Fig. 14 right panel).

BepiColombo and other missions

At present time, the semi-major axis of the two orbiters (MPO & MMO) of the Bepi-Colombo mission are respectively equal to 3 394 km and 8 552 km. The MPO altitude corresponds to our (E') phase space where the eccentricity vector has a circular concentric motion. The MMO initial conditions, without thrust correction, leads to a crash onto the Mercury surface after 3 years. Based on our theory (and the underlying assumptions), we can choose another initial condition a = 7 355 km and e = 0.652, that avoids forever the crash on Mercury and whose eccentricity vector is fixed.

Conclusions

The orbit dynamics of a space probe orbiting a planet or a natural planetary satellite has been investigated. The proposed model includes the effects of J 2 for the central body and the perturbation of the third body. We have developed a doubly averaged Hamiltonian and studied the location of the stable and unstable frozen orbits. Our method allows to have an analytical global approach to search all these equilibria and to compute the periods of the free librations at the equilibria. The analytical results have been checked and validated numerically by performing numerical integrations of the complete systems. Our theory is able to explain the behavior of our preliminary numerical investigations where the variation of the amplitude of the eccentricity is null and a separatrix is detected. The theory is general enough to be applied to a wide range of probes around any planet or any natural planetary satellite, provided that they respect the hypotheses used to obtain our Hamiltonian model.

We have shown the protection mechanism of the coefficient J 2 on the increasing of the eccentricity due to Kozai-Lidov effect. This mechanism is therefore able to produce a larger number of frozen orbits than those found in the only Kozai-Lidov problem. We have also explained the asymmetry of the frozen equilibria caused by the addition of the coefficient J 3 . We have also brought to the light a local deformation of the Kozai-Lidov equilibria that appears at high eccentricity, high inclination and large value of γ.

These equilibria and their periods can be modified by including higher order terms of gravity field. However they will not change drastically the behavior of our phase space. Henceforth, the frozen orbits presented here provide good initial conditions to find periodic or frozen orbits in a osculating full model.

This theory provides analytical results to justify the choice of the intial semi-major axis and eccentricity of an orbiter for future missions around planets or planetray satellites.

Interested (separated) improvements would be to include the C 22 coefficient of the central body or the inclination of the third body. In these cases the ascending node of the probe would take place in the Hamiltonian leading to an extra degree of freedom i.e. 4D phase plots.

  Fig.1The eccentricity computed as a function of the initial eccentricity e 0 and the initial semi-major axis a 0 . The equations of motion include the central body attraction, the harmonics J 2 , C 22 , J 3 , the solar interaction as well as the perturbing effects of the solar radiation pressure (A/m = 0.01m 2 /kg). The eccentricity step is 0.005 and the semi-major axis step is 35 km. The initial conditions are i 0 = 90 • , Ω 0 = 67.7 • , ω 0 = -2 • and M 0 = 36.4 • . The integration time is 200 years from epoch fixed at 14 September 2019. The patterns have been obtained by plotting the amplitude of variation of the eccentricity in left panel and the logarithm of the second derivative, log | ∂ 2 ν ∂k ∂h |, of the fundamental frequency of the eccentricity vector evolution in right panel.

Fig. 2

 2 Fig.2Reference planes, for the general theory on the left and for our particular case on the right.

Fig. 3

 3 Fig.3Relation between γ and the semi-major axis of a test particule orbiting the central body (terrestrial planets or Europa). The third body are respectively the Sun and Jupiter.

Fig. 4

 4 Fig. 4 Values of G at Kozai-Lidov stable (Eqs. 18 and 19) equilibria (vertical equilibria: ω = ±π/2) computed as a function of H 2 and γ. These equilibria are always stable. The color code indicates the value of G at the equilibrium.

Fig. 5

 5 Fig.5Value of G at stable (on the left) (Eqs. 24 and condition 26) equilibrium and unstable (on the right) (Eqs. 24 and condition 27) equilibrium, computed as a function of H 2 and γ. The color code indicates the value of G at the equilibrium (horizontal equilibria: ω = 0, π). The inset on the left panel shows the same plot of left panel but using a wider color code.

Fig. 7

 7 Fig.7Examples of some generic contour plots of the Hamiltonian (13) in (k, h, ı) space for each region of Fig.6. The inclination ı is given in degrees and the semi-equinoctial elements (k, h) are given by (k,h) = ( √ 1 -G 2 cos ω, √ 1 -G 2 sin ω).The green point and red cross are respectively the stable and unstable points. In the polar projection, the radius is the eccentricity e and the angle is the pericenter ω in degrees.

Fig. 9

 9 Fig.9Plot of the periods at the stable 19) and (0, 0) equilibrium (with conditions 31) respectively in the upper, center and lower panels. The color code indicates the period (truncated to 100 years) of the fundamental frequency at the equilibrium. On the left panels, the period with respect to the parameters (γ, H 2 ). On the right panels, the location of the stable equilibrium in the physical space (a, e, ı) with its period. For the equilibrium (0, 0), e is always equal to 0.

Fig. 10

 10 Fig.10Comparison between analytical and numerical results. For the left and right panel, the initial conditions are a 0 = 6 407 km (γ = 1.000296), Ω 0 = M 0 = 0 • . In the left panel, for the lower orbit, we take e 0 = 0.545055 (G = 0.8384), ı 0 = 76.646989 • (H 2 = 3.7492 × 10 -2 ) and ω 0 = 180 • ; for the upper orbit, we take e 0 = 0.6 (G = 0.8), ı 0 = 78.221768 • (H 2 = 0.26666666) and ω 0 = 0 • . For the right panel, the initial conditions are e 0 = 0.01 (G = 0.99995), ı 0 = 69.73104 • (H 2 = 0.1199999) and ω 0 = 0 • . For the middle panels, the initial conditions are, for the lower panel a 0 = 4 650 km (γ = 0.17096), e 0 = 0.3 (G = 0.9539392), ı 0 = 77.89775 • (H 2 = 0.04) and Ω 0 = M 0 = ω 0 = 0 • ; for the upper panel, we take the MMO initial conditions: a 0 = 3 394 km (γ = 0.0417), e 0 = 0.163229 (G = 0.986), ı 0 = 90 • (H 2 = 0) and Ω 0 = 67.7 • , ω 0 = 16 • , M 0 = 0 • . The numerical model takes into account the contribution of J 2 and 22 and the solar gravitational effect, with starting epoch fixed at 14 September 2019. The analytical model is based onEquations (14, 16, 17). We plot the numerical integrations with continued lines and the analytical results with dashed lines. The time of the integrations is equal to 180 years. In the right panel, the numerical integration leads to a crash onto the planet after 40 years.

Fig. 11

 11 Fig. 11 Comparison between analytical and numerical results for exact polar orbiter. The initial conditions are a 0= 6 000 km (γ ≃ 0.72), ı 0 = 90 • , Ω 0 = 67.7 • , ω 0 = -2 • , M = 36.4 • .The numerical and analytical model are the same of Figure10. In dashed line the analytical result and in continued line the numerical integration. On the right a blow-up of the center of libration.

Fig. 13

 13 Fig. 13 Local deformation of the Kozail-Lidov equilibrium. The bifurcation lines in the left panel with the new region (G) demarcated by the two curves KL 3 . Example of generic contour (for (G) region) of the Hamiltonian (13) in (k, h, ı) space in the middle panels. A zoom of the local deformation in the right panels.

7. 3 . 1

 31 Vertical equilibria -Kozai-Lidov equilibria: cos ω = 0 ⇔ ω = ±π/2

Fig. 14

 14 Fig.14Distortion of the phase space (for a Mercury's orbiter) due to J 3 effect. The initial conditions are a = 5 900 km (γ ≃ 0.66) and H 2 = 0.06. J 3 (respectively δ) is equal to 0 (0), -J 2 /10 (-0.041356) and -J 2 /2 (-0.20678) in left, center and right panels.

  16 γ 6 + 2 963 520 H 12 -1 024 H 10 γ 4

	3087 and H 2 ≤ 1 + 1 512 630 H 8 -13 965 H 6 -22 235 661 H 10 γ 2 + 12 = 0	(20

Table 2

 2 Comparison between the period of the equilibria determined in the analytical model and the period numerically NAFF. The Horiz. ⋆ corresponds to the equilibrium in the blow-up of Fig.11

		Initial condition		Period [year]	Error
	Which	a	e	ı	Analytical Numerical relative
	equi.	[km]		[degree]			%
	Kozai	5 750 0.4731	58.37	29.30	29.25	0.17
	Horiz.	8 083 0.4922	77.68	35.67	35.61	0.17
	Horiz.	5 818 0.5418	71.93	42.17	42.26	0.21
	Horiz. ⋆ 6 000	0.369	90.00	44.576	44.528	0.11
	(0,0)	3 429	0.0	47.64	9.127	9.135	0.08
	(0,0)	4 731	0.0	77.01	56.594	55.274	2.38

  000 H 16 γ 6 + 2 963 520 H 12 -1 024 H 10 γ 4

	+ 1 512 630 H 8 -13 965 H 6 -22 235 661 H 10 γ 2 + 12 = 0 and H 2 ≤ 1 3087

This can be done calculating the discriminant of the polynomial given by the first eq. 18 to find the roots (G).

This region is the (G) zone in Figure13.

At this inclination, due to J 2 effect, the argument of perigee remains nearly constant for a long period of time. Molniya orbits are named after a series of Soviet/Russian Molniya communications satellites which have been using this type of orbit since the mid 1960s.

Acknowledgements The authors thank B. Noyelles and A. Albouy for fruitful discussions, the IMCCE team for their welcome and B. Meyssignac for initializing discussion.

Numerical simulations were made on the local computing ressources (Cluster URBM-SYSDYN) at the University of Namur (FUNDP, Belgium).

This work was partly supported by the fellowship Concours des bourses de voyage de la Communauté Française de Belgique obtained by Nicolas Delsate.

Finally, the authors thank the two referees for their suggestions that contributed to improve this paper.


Nicolas Blondiaux
email: nicolas.blondiaux@chru-lille.fr

Frédéric Wallet

Raphaël Favory

Thierry Onimus

Saad Nseir

René J Courcol

Alain Durocher

Micheline Roussel-Delvallez

Daily serum piperacillin monitoring is advisable in critically ill patients

Keywords:

A c c e p t e d M a n u s c r i p t

ABSTRACT

The aim of the present study was to evaluate the benefit of monitoring serum piperacillin concentrations in critically ill patients. This was an 11-month, prospective, observational study in a 30-bed Intensive Care Unit in a teaching hospital, involving 24 critically ill patients with evidence of bacterial sepsis. All patients received a 66 mg/kg intravenous bolus of piperacillin in combination with tazobactam (ratio 1:0.125) followed by continuous infusion of 200 mg/kg/24 h. The dosage was adjusted when the serum piperacillin concentration either fell below 4 the drug's minimum inhibitory concentration (MIC) for the causative agent or exceeded the toxic threshold of 150 mg/L. With the initial regimen, serum piperacillin concentrations were within the therapeutic target range in only 50.0% of patients (n = 12). This proportion increased to 75.0% (18 patients) (P = 0.006) following dosage adjustment. For patients with low initial serum piperacillin concentrations (n = 8), the percentage of time during which the concentration remained above 4 MIC (%T>4 MIC) was 7.1 ± 5.9% before dosage adjustment and 27.3 ± 8.6% afterwards. In conclusion, in critically ill patients, monitoring and adjustment of serum piperacillin levels is required to prevent overdosing and might also help to correct underdosing, an important cause of antibiotic therapy failure.

Introduction

Broad-spectrum -lactams are commonly used to treat severe infections (either alone or combined with other antibiotics), especially in the Intensive Care Unit (ICU) setting. The bactericidal effects of these antibiotics are mostly time-dependent. Laboratory and clinical evidence indicates that to achieve effective bactericidal levels, a serum concentration of at least 4 times the minimum inhibitory concentration (MIC) for the pathogen involved may be required [START_REF] Craig | Continuous infusion of -lactam antibiotics[END_REF][START_REF] Macgowan | Pharmacodynamics of antimicrobial agents and rationale for their dosing[END_REF]. A related determinant of -lactam efficacy is the time during which the drug's serum concentration remains above the MIC threshold (T>MIC) [START_REF] Craig | Pharmacokinetic/pharmacodynamic parameters: rationale for antibacterial dosing of mice and men[END_REF]. Continuous infusion of these antimicrobial agents facilitates the optimisation of their pharmacodynamic profiles [START_REF] Kasiakou | Continuous versus intermittent intravenous administration of antibiotics: a meta-analysis of randomised controlled trials[END_REF][START_REF] Roberts | Continuous infusion of lactam antibiotics in severe infections: a review of its role[END_REF].

In ICU patients, it is more difficult to maintain stable serum antibiotic levels. Apart from the excretory impairments observed in these patients, several other phenomena can significantly alter drug distribution and/or elimination profiles [START_REF] Pea | Antimicrobial therapy in critically ill patients: a review of pathophysiological conditions responsible for altered disposition and pharmacokinetic variability[END_REF]. Furthermore, ICU patients show great interindividual variability in their pharmacokinetic parameters [START_REF] Mehrotra | Antibiotic pharmacokinetic and pharmacodynamic considerations in critical illness[END_REF], which can also change quickly in response to improvements in or worsening of the clinical condition. Monitoring of serum antibiotic concentrations and subsequent dosage adjustment can avoid underdosing or overdosing and limit the risk of therapeutic failure and toxicity.

Piperacillin/tazobactam (PIP/TAZ) is a -lactam/-lactamase inhibitor combination with activity both against Gram-negative and Gram-positive bacteria (including -lactamase Over the last decade, new modes of administration have significantly optimised the antimicrobial activity of PIP/TAZ. In several pharmacokinetic studies, continuous infusion of PIP/TAZ was superior to intermittent injection (despite the use of lower doses) [START_REF] Buck | Pharmacokinetics of piperacillin-tazobactam: intermittent dosing versus continuous infusion[END_REF]9] and also appeared to be a better pharmacoeconomic option [START_REF] Kotapati | The clinical and economic benefits of administering piperacillin-tazobactam by continuous infusion[END_REF]. Under these supposedly optimal pharmacokinetic conditions, the objective of this study was to assess the utility of daily serum piperacillin monitoring to achieve pharmacodynamic targets associated with optimal piperacillin activity.

Patients and methods

Subjects

Study subjects were adult patients admitted to the Department of Respiratory Emergencies and Intensive Care at Lille University Medical Center (Lille, France) with various pathologies (Table 1). The study's main inclusion criterion was the presence of a community-or hospital-acquired infection requiring antibiotic therapy with PIP/TAZ.

Patients with known hypersensitivity/intolerance to penicillins or a previous infection with PIP/TAZ-resistant bacteria were excluded from the study.

A c c e p t e d M a n u s c r i p t 6

Antibiotic administration

All patients received PIP/TAZ via continuous infusion. The regimen consisted of a loading bolus of 66 mg/kg piperacillin and 8.25 mg/kg tazobactam (Tazocilline ® ; Wyeth Pharmaceuticals, Paris, France) injected over 30 min followed by constant-rate infusion of 200 mg/kg/24 h piperacillin and 25 mg/kg/24 h tazobactam. Although this dosage may not be adapted to patients weighing more than 120 kg, we did not use daily doses above 24/3 g of PIP/TAZ because they are not recommended by the manufacturer. Furthermore, none of the enrolled patients weighed more than 105 kg. In patients with impaired renal function [creatinine clearance (CL Cr ) ≤ 30 mL/min, calculated according to the Cockcroft and Gault formula], the maintenance dose was halved.

Target serum concentrations

Dosage adjustment was performed when the serum piperacillin concentration was below 4 MIC for the identified pathogen or >150 mg/L (defined as the toxicity threshold, based on empirical data). Above this concentration, the likelihood of neurological adverse events increased (personal data). When the pathogen was not known, the lower target value was defined as 4 the critical concentration above which all bacterial strains are considered to be resistant to piperacillin, i.e. 4  16 mg/L = 64 mg/L (according to guidelines issued by the French Microbiology Society's Antibiotics Committee).

A c c e p t e d M a n u s c r i p t 7

Analytical method

Serum piperacillin concentrations were determined by high-pressure liquid chromatography using a 322-A pump, an autoinjector, a 332 variable ultraviolet wavelength detector (all from Serlabo Technologies, Entraigues-sur-la-Sorgue, France) and a C18 Resolve column (4 m, 3.9  150 mm; Waters Associates, Milford, MA). The mobile phase consisted of 0.02 M ammonium acetate buffer and acetonitrile in a 76:24 ratio (v/v) at pH 5.0. The flow rate was 1 mL/min. Pooled human serum was used to prepare standards, to check samples and to dilute serum samples, as required. The assay was linear over the range 1 mg/L to 200 mg/L. Only the free fraction of piperacillin was measured, as an initial step of protein precipitation with acetonitrile was carried out [START_REF] Rudrik | Determination of penicillinase-resistant penicillins in serum using high pressure liquid chromatography[END_REF].

Blood sampling

Blood samples for pharmacokinetic analysis were obtained from a peripheral venous catheter placed in the arm not used for drug infusion. Blood was collected before starting drug administration (t = 0) and 1.5 h afterwards (t = 1.5) to ensure that the bolus injection had been completed. A sample was then taken each morning at steady state until the end of the course of treatment.

Adverse events

Convulsions, changes in white blood cell count, modification of CL Cr and signs of cutaneous intolerance were monitored throughout the study.

Statistical analysis

A paired Student's t-test was used to check for statistical significance. The serum piperacillin concentration was monitored in 4 to 21 samples per patient (mean  S.D., 10 ± 4.6). With the initial regimen, the serum piperacillin concentrations were above the previously defined target value (i.e. >4 MIC) in 12 (50.0%) of the patients (Patients 1-6, 9, C, D, G, H and I) (Table 1). This goal was achieved in six additional patients (Patients 7, 8, 10, A, E and F) following dosage adjustment (P = 0.006). In four patients (Patients 8, A, E and F), antibiotic levels were considered toxic A c c e p t e d M a n u s c r i p t 9 (>150 mg/L), thus the dosage was reduced and no adverse events were observed.

Results and discussion

Conversely, a dosage increase in two patients (Patients 7 and 10) of the eight patients with below-target serum piperacillin concentrations enabled achievement of the therapeutic value. Dosage adjustment was not efficient in the six remaining patients (Patients 11-15 and B). Despite many dosage changes (with up to 18 g of piperacillin per day in some cases), serum piperacillin concentrations >4 MIC were never achieved. Interestingly, all but one (Patient B) of these subjects had extensive cellulitis, which might have influenced the pharmacokinetics of the antibiotic [START_REF] Gould | Piperacillin/tazobactam in the treatment of serious acute soft tissue infection[END_REF]. In this clinical setting, subcutaneous inflammation might generate an additional distribution compartment and could be responsible for extracellular and extraplasmatic leakage of the antibiotic. This hypothesis is supported by the substantial increase in the volume of distribution of PIP/TAZ seen in burns victims compared with healthy individuals (indicating translesional diffusion of the antibiotic) [START_REF] Bourget | Clinical pharmacokinetics of piperacillin-tazobactam combination in patients with major burns and signs of infection[END_REF]. Further investigation is now required to establish suitable antibiotic regimens in this setting; nevertheless, daily monitoring of serum piperacillin concentrations significantly increased achievement of the therapeutic target concentration.

By contrast to that shown for fluoroquinolones [START_REF] Drlica | Mutant selection window hypothesis updated[END_REF], the possibility of selecting resistant mutants at serum piperacillin concentrations <4 MIC is not supported by experimental data. However, several authors have suggested that prolonged percentage of time during which the concentration remains above 4 MIC (T>4 MIC), instead of T>MIC, is an important pharmacokinetic/pharmacodynamic (PK/PD) parameter that better predicts treatment outcome with -lactams (especially for critically ill patients) [START_REF] Macgowan | Pharmacodynamics of antimicrobial agents and rationale for their dosing[END_REF][START_REF] Roberts | Continuous infusion of lactam antibiotics in severe infections: a review of its role[END_REF][START_REF] Mouton | Killing of Pseudomonas aeruginosa during continuous and intermittent infusion of ceftazidime in an in vitro pharmacokinetic model[END_REF]. Here, before antibiotic dosage adjustment and 27.3 ± 8.6% afterwards (Fig. 1). This significant difference (P = 0.03) demonstrates the beneficial impact of regular serum piperacillin monitoring on these important therapeutic PK/PD parameters.

Conclusion

In critically ill patients, monitoring and adjustment of serum piperacillin levels are warranted to prevent overdosing and to help correct underdosing, an important cause of antibiotic therapy failure. Additional trials in a more homogeneous patient population are required to underpin these preliminary findings.

A c c e p t e d M a n u s c r i p t Characteristics of the 24 patients treated with piperacillin/tazobactam (PIP/TAZ) (a) Characteristics of the 15 patients with non-documented infection (no pathogen isolated). The lower target value was defined as 4 the critical concentration above which all bacterial species are considered to be resistant to PIP, i.e. 4  16 mg/L = 64 mg/L (according to guidelines issued by the French Microbiology Society's Antibiotics Committee). For these patients, target values are between 64 mg/L and 150 mg/L (toxicity threshold); otherwise dosage was changed Patien t no.

  The broad spectrum of action of PIP/TAZ makes it an appropriate choice in the ICU.

  ) and included Pseudomonas aeruginosa (Patients A and B), Escherichia coli (Patients C and D), Klebsiella pneumoniae (Patients E and F), Acinetobacter baumannii (Patient G), meticillin-susceptible Staphylococcus aureus (Patient H) and Enterobacter aerogenes (Patient I). In these strains, the MIC of piperacillin (when combined with tazobactam) ranged from 4 mg/L to 16 mg/L (mean value 7.6 ± 5.1 mg/L).

  were monitored daily for between 3 days and 20 days (mean  S.D., 9 ± 5.2 days). For the eight patients (Patients 7, 10-15 and B) with low initial serum piperacillin concentrations, the mean percentage T>4 MIC was 7.1 ± 5.9%

Fig. 1 .

 1 Fig. 1. Impact of piperacillin/tazobactam dosage changes on

( b )

 b Characteristics of the 9 patients with documented infection. The lower target value was defined as 4 the PIP MIC for the identified pathogen. For these patients target values are between 4 the PIP MIC and 150 mg/

u s c r i p t 5 HAP

 5 , hospital-acquired pneumonia; sUTI, severe urinary tract infection; P.a., Pseudomonas aeruginosa; E.c., Escherichia coli; K.p., Klebsiella pneumoniae; A.b., Acinetobacter baumannii; S.a., Staphylococcus aureus; E.a., Enterobacter aerogenes. a Creatinine clearance 30 mL/min, calculated according to the Cockcroft and Gault formula. b PIP/TAZ combination. c Serum PIP concentration (free fraction): mean ± S.D. of 2 to 17 samples according to the patient.

Table 1b

 1b In total, 24 patients [22 males and 2 females; mean ± standard deviation (S.D.) age, 62.5 ± 14.3 years; mean ± S.D. weight, 78.7 ± 16.8 kg] were initially included. Patient characteristics and PIP/TAZ dosage management are presented in Table 1. Seven patients (2, 3, 8,12, A, E and F; Table 1) displayed severe renal function impairments (CL Cr ≤ 30 mL/min) The infecting pathogen was only documented in 9 of the 24 cases (

Funding

None.

Competing interests

None declared.

Ethical approval

This study was conducted in accordance with the Declaration of Helsinki, and written informed consent was obtained from all patients or their legal guardians. Edited Figure 1